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1 Introduction 

As a tool for large-scale computations in fluid dynamics, spectral methods were prophesied in 
1944, born in 1954, virtually buried in the mid-1960’s, resurrected in 1969, evangalized in the 
1970’s, and catholicized in the 1980’s. The use of spectral methods for meteorological prob- 
lems was proposed by Blinova in 1944 and the first numerical computations were conducted 
by Silberman (1954). By the early 1960’s computers had achieved sufficient power to permit 
calculations with hundreds of degrees of freedom. For problems of this size the traditional way 
of computing the nonlinear terms in spectral methods was outrageously expensive compared 
with finite-difference methods. Consequently, spectral methods fell out of favor. The expense 
of computing nonlinear terms remained a severe drawback until Orszag (1969) and Eliasen, 
Machenauer, and Rasmussen (1970) developed the transform methods that still form the back- 
bone of many large-scale spectral computations. The original proselytes of spectral methods were 
meteorologists involved in global weather modeling and fluid dynamicists investigating isotropic 
turbulence. The converts who were inspired by the successes of these pioneers remained, for the 
most part, confined to these and closely related fields throughout the 1970’s. During that decade 
spectral methods appeared to be well-suited only for problems governed by ordinary differential 
equations or by partial differential equations with periodic boundary conditions. And, of course, 
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the solution itself needed to be smooth. Some of the obstacles to wider application of spectral 
methods were: 1) poor resolution of discontinuous solutions; 2) inefficient implementation of 
implicit methods; and 3) drastic geometric constraints. All of these barriers have undergone 
some erosion during the 1980’s, particularly the latter two. As a result, the applicability and 
appeal of spectral methods for computational fluid dynamics has broadened considerably. 

The motivation for the use of spectral methods in numerical calculations stems from the 
attractive approximation properties of orthogonal polynomial expansions. Suppose, for example, 
that a function, ti(x), is expanded in a Legendre series on [-1,1]: 

«*(*) = UnM 1 ). (!) 

n~0 

where L n (x) is the Legendre polynomial of degree n. The classical form of the expansion 
coefficients (or spectra) is 

i r 1 

«n = (*+-) J i u(x)L n (x)dx, for n > 0. (2) 

An integration-by-parts argument reveals that 

n p ii n — ► 0 as n — > oo, for all p> 0, (3) 

provided that u is infinitely differentiable. Consequently, the approximation error decreases 
faster than algebraically. This rapid convergence is referred to as infinite-order accuracy, expo- 
nential convergence, or spectral accuracy. Our primary concern in this review is on numerical 
methods for partial differential equations that exhibit spectral accuracy for infinitely differen- 
tiable solutions. 

The two most common types of finite approximations are based on truncation and interpo- 
lation. The Legendre series truncated after N terms is 

AT 

P N u(x) = J2 u n L n {x). (4) 

n= 0 

It forms the foundation of Legendre Galerkin and tau methods. The Legendre interpolant of 
u(x), on the other hand, is used for collocation approximations, and is given by 

N 

I N u(x) = ^2 u n L n {x), (5) 

n= 0 

where the coefficients, u n , of the interpolant are determined by the condition, 


I N u(xj) = u(xj), for j = 0, 1, . . . , N, 


( 6 ) 


special, so-called collocation, points in [-1,1], denoted by Xj. For most problems, the optimal 
choice of these collocation points is 


x 


j 


' -1 j = 0 

< j'th zero of L' N j = 1, 2, . . . , IV — 1 . 
k +1 J = N 


( 7 ) 


This choice of collocation points yields an extremely accurate approximation to the integral 
appearing in Eq. (2): 


1 * 
u n = — 


7n j=0 


for j = 0,1,..., N, 


( 8 ) 
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where the Gauss-Lobatto weights, 


2 1 

“ v - mir + i)LU*i)’ 

and the normalization constants, 


for j = 0,1 


In = 


( 



0<n< N - 1 
n = N 


( 9 ) 


( 10 ) 


Whether (2) or (6) is used for the expansion coefficients, the expansion (l) is differentiated 
analytically to form the approximations to whatever derivatives are required for the problem at 
hand. 

A graphical distinction between traditional approximations and spectral ones is provided in 
Figure 1 for the simple task of estimating the derivative of the function 1 + 8tn(2nx + tt/ 4) on 
[-1,1] from the values of the function at a finite number of grid-points. A finite-difference or 
finite-element approximation uses local information to estimate derivatives, whereas a spectral 
approximation uses global information. In this figure a second-order (central) finite- difference 
approximation is compared with a Legendre spectral collocation approximation. The finite- 
difference approximation estimates the derivative at, say, x = 0, from the parabola which 
interpolates the function at x = 0 and the two adjacent grid-points. A separate parabola is 
used at each grid-point. The spectral approximation, on the other hand, uses all the available 
information about the function. If there are N + 1 grid-points, then the interpolating polyno- 
mial from which the derivative is extracted has degree 1NT, and the same polynomial is used for 
all the grid-points. For small N the accuracy of the spectral and the finite-difference results 
are comparable. However, as N increases, the accuracy of the spectral approximation increases 
dramatically (its error decays exponentially), whereas the finite-difference result improves only 
algebraically. 

The high accuracy achieved by the spectral method for a modest number of degrees of free- 
dom on this simple one-dimensional approximation problem has been demonstrated repeatedly 
for multi-dimensional, nonlinear problems in fluid dynamics. This is particulary true for those 
computations, such as the numerical simulation of transition to turbulence, for which great accu- 
racy is essential rather than a luxury. Figure 2, taken from a spectral simulation by Streett and 
Hussaini (1989), exemplifies this. The cascade of bifurcations leading to an aperiodic, weakly 
turbulent (chaotic) state in Taylor-Couette flow is an established experimental fact, but has 
proven to be an elusive goal for numerical simulations. In order to capture properly the dynam- 
ics of such states, long-time integrations of the Navier-Stokes equations are required. Numerical 
errors which are seemingly innocuous over a single time-step or in a steady-state result, never- 
theless can accumulate to dangerous levels over long times. Hence, highly accurate calculations 
are imperative. The aperiodic time history and the resultant broadband frequency spectrum 
shown in the figure reflect the high accuracy achieved by the spectral simulation. 

One of the objectives of these notes is to provide a basic introduction to spectral methods 
with a particular emphasis on applications to computational fluid dynamics. Another objective 
is to summarize some of the most important developments in spectral methods in the last two 
years. The fundamentals of spectral methods for simple problems will be covered in depth, and 
the essential elements of several fluid dynamical applications will be sketched. The recent book 
by Canuto, Hussaini, Quarteroni, and Zang (1988) (referred to hereafter as CHQZ) contains 
a comprehensive discussion of both the implementation and the rigorous numerical analysis of 
spectral methods. The first part of these notes will rely heavily upon that text, with frequent 
references to it for details, and will also include some material from our earlier reviews (Hussaini, 
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Streett, and Zang 1983, Hussaini and Zang 1987). Moreover, the notation in these notes conforms 
closely to that used in the text. Since CHQZ contains an exhaustive list of references through 
1986, we shall not bother here to cite the original sources for most of the methods that we discuss. 
Instead, the bibliography for these notes focuses on those references which have appeared since 
1986 and which pertain to new developments in spectral methods rather than just additional 
applications of established algorithms. 

Additional general references on spectral methods are the monographs by Gottlieb and 
Orszag (1977) and Mercier (1981), the proceedings of the ICASE Spectral Methods Workshop 
edited by Voigt, Gottlieb, and Hussaini (1984), and doubtless also the forthcoming proceedings 
of the International Conference on Spectral and High-Order Methods for Partial Differential 
Equations to be held in Como, Italy in June, 1989. 


2 Spectral Approximations 

This chapter provides a basic introduction to Fourier and Chebyshev spectral approximations. 
Fourier approximations to infinitely differentiable, periodic functions exhibit exponential con- 
vergence. They do not behave nearly so well for non-periodic functions, regardless of their 
differentiability. For such functions expansions in Jacobi polynomials are appropriate. Cheby- 
shev polynomials are the most useful polynomials in this class. (Legendre polynomials also 
belong to this class, but they suffer from the lack of a fast transform.) 

The focus of this chapter is on practical details. It closely parallels the discussion in CHQZ, 
Ch. 2. The extensive literature on the rigorous numerical analysis of spectral approximations 
is surveyed in CHQZ, Ch. 9. 


2.1 Fourier Expansions 

The continuous Fourier expansion of a function, u(x), on (0, 2x) is 

«(*)= £ 

*=-00 


( 11 ) 


( 12 ) 

(13) 


where the Fourier coefficients, 

The orthogonality condition is 

1 f 2 * ■ 

± e* px dx = 8po, 

2n Jo 

where S^i denotes the Kronecker delta-function. 

Let us begin by analyzing the convergence of this series for a function which is infinitely 
differentiable and which, together with all its derivatives, has period 2ic . Using integration-by- 
parts we obtain 

<“ = - sse - m *— * + hrJo'^‘"" dx - (14) 

Invoking periodicity, 

(15) 


The Riemann-Lebesgue lemma implies that 

ii k = 0(k ~ 2 ) 


as k — > oo. 


(16) 


5 


Repeated integrations-by-parts leads to 




{ik) 

hence to the conclusion that, for all positive integers p, 

ilk = 0(k~ p ~ 1 ) as k — ♦ oo. 

Thus, 


k p Uk — » 0 

for all positive integers p. 

The truncated Fourier series is defined to be 


as k —f oo 


(17) 

(18) 
(19) 


N 


-1 


P N u(x) = ^2 u k e' kz . 


( 20 ) 


k=-H 

2 


The property (19) implies that the truncated Fourier series converges faster than any finite 
power of 1 fN. In theoretical discussions of Fourier series, the truncation is customarily defined 
by 

N_ 

Pnu(x) = £ u k e ik *, (21) 

rather than by (20). However, because of the properties of the Fast Fourier Transform (see 
below), the truncation (20) is the one that is invariably used in practice. 

For a less well-behaved function the repeated integrations-by-parts must stop at some point, 
either because u(x) has only finitely many derivatives or because for the q th derivative uM(27r“) ^ 
u^(0 + ). In such a case the Fourier coefficients, and hence also the truncated Fourier series, 
decay only algebraically. 

The discrete Fourier expansion is closely related. It is inherently a finite series, and is 
typically based on values of t*(x), denoted by uy, which are given at the collocation points, 

Xj = 2xj/ N , J — 0,1,...,# — 1. (22) 


The discrete Fourier coefficients are 


1 jV ‘ 1 

^=nT, 

j = 0 





The inversion formula is 


and the orthogonality condition is 





k= 


N 

2 


J_ y^ 1 e i P x, i= f 1 p=Nm,m = 0,±1,±2,... 

N " ) 0 otherwise 

j=o t 


(23) 


(24) 


(25) 
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Associated with the discrete Fourier series is the interpolating trigonometric polynomial, 


fr-i 


/W*)= E 


L N 


which is defined for all a: € [0, 2*r]. Note that 

Uj = u(xy) = /ivw(xy). 


(26) 


(27) 


It is customary to refer to u(z) and uy as the physical space representations of u, and to u* 
and u k as its transform space representations. 

For a discrete Fourier series, both the function itself and its Fourier coefficients are periodic: 

Uj+N = «y ( 28 ) 

U/fc+AT = ( 29 ) 

For a continuous Fourier series, of course, only the function itself is periodic. 

The sums in (23) and (24) can be accomplished by means of the Fast Fourier Transform 
(FFT). The simplest FFT requires N to be a power of 2. If the data are fully complex it 
requires 5 N log^N — 6 N real operations, where addition and multiplication are counted as 
separate operations. In most applications, u is real, and therefore u* = u*_ k . In this case 
the operation count is halved. Fast Fourier Transforms which allow factors of 2, 3, 4, 5, and 
6 are widely available (Temperton 1983), and they offer a 10-20% reduction in the operation 
count over the basic power of 2 FFT. For simplicity, we shall often use just 5 N log^N as 
the operation count for a complex FFT. Thus, to compute the discrete Fourier coefficients of 
Uj requires log^N operations; the cost for computing the grid-point values, tiy, given the 
discrete Fourier coefficients is the same. A more complete discussion of FFT’s (including several 
FORTRAN programs) is available in CHQZ, Appendix B. 

The FFT’s that are typically used in spectral methods require N to be even. This accounts 
for the use of (20) rather than (21) for the truncated Fourier series. The y mode is something 
of a nuisance (see also the remark after (37) below), and it is probably advisable in practice 
simply to leave it out of the finite Fourier representation altogether. 

The discrete Fourier coefficients are related to the continuous ones by 

oo 

Ufc = U* + “k+Nrn 


k = 


N N 


2 ’ 2 


+ 1 , . • 


N 


- 1. 


(30) 


m^0 

This expression shows that the k th mode (or Fourier coefficient) of the trigonometric interpolant 
of u depends not only on the k th mode of u, but also on all the modes of u which “alias” the k th 
one on the discrete grid. The (k + Nm) th frequency aliases the k th frequency on a grid consisting 
of N points; they are indistinguisable at the nodes since e , ( k + Nm ) x i — e ikx > . This phenomenon 
is illustrated for N = 8 in Figure 3. Shown there are three sine waves with frequencies k = 6, 
-2, and -10. Superimposed upon each wave are the eight grid-point values of the function. In 
each case these grid-point values coincide with the k = -2 wave. 

An equivalent formulation of (30) is 

Inu = Pnu + Rnu, (31) 
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with 


n_i 

R N“= E ( E Wm)« ita . (32) 

jfe— — — 

*— ~2 m=- oo 

m^O 

The error, ffo u, between the interpolating polynomial and the truncated Fourier series is called 
the “aliasing error” . It is orthogonal to the truncation error, tt — P/yu, so that 

II u - I N u || 2 =|| u - P N u || 2 + || R n u || 2 . (33) 


Hence, the error due to interpolation is at least as great as the error due to truncation. 


In a Fourier spectral method the derivative of a function, u(z), is approximated by the 
analytic derivative of the finite series which approximates u(z). In the case of the truncated 
Fourier series, 


djPyu) 

dx 


f- 1 

E »W fcz . 


(34) 


This is called the Fourier Galerkin derivative. Note that 


d{P N u) 

dx 




(35) 


i.e., truncation and differentiation commute. 
Likewise, for the interpolating Fourier series, 


dx 


f- 1 

E *u k e ik *. 


(36) 


This is called the Fourier collocation derivative. In this case, however, 


dx T 


«=)• 


(37) 


i.e., truncation and differentiation do not commute. Incidentally, the y term in the sums (34) 
and (36) make a purely imaginary contribution to the derivative and are best left out entirely. 

In a typical spectral Galerkin method, the fundamental representation is in transform space. 
Differentiation is accomplished by a simple multiplication. For the spectral collocation method, 
on the other hand, the fundamental representation is in physical space. In this case, Fourier 
differentiation using a transform method (employing two FFT’s) is accomplished in three steps: 
(1) compute the discrete Fourier coefficients via (23); (2) compute the derivative coefficients 
fiW = tkujc] (3) perform the sum (24) with uj^ in place of u*. Steps (1) and (3) require 
| N log 2 N — 3 N operations each, and step (2) costs N multiplications. Thus, the total cost of 
the derivative is 5 N log 2 N - 5 N operations. 

Fourier collocation differentiation can be represented by a matrix: 


djlNu) 

dx 


(xjt) = D$u h 


1=0 


where 


£)(J) _ / |(-i)* + M ! V r ) 


l^k 

l = k ' 


(38) 


(39) 
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For the second derivative we have 


where 


<P(Itfu) f ^ n( 2 ) 

dx 2 (**) = 2^ D kl«l> 
aX 1=0 

(40) 

\ -(iV 2 -2)/12 / = k 

(41) 


Fourier collocation derivatives obviously may be calculated by simply performing the matrix- 
multiplies implied by (38) or (40). The operation count for this approach is 2A r2 . For large N, 
this is much larger than the operation count for the transform method. Moreover, the matrix- 
multiply approach is more susceptible to round-off errors. For small or moderate values of N, 
however, the matrix-multiply approach is competitive. Some timings are available in CHQZ, 
Sec. 2.1.3. 


2.2 Chebyshev Expansions 

Spectral methods for non-periodic problems are typically based upon Chebyshev polynomials. 
These are defined on [-1,1] by 

T k (x) = cos(A; co8~ 1 x ) A; > 0. (42) 

The first several are To(x) = l,Ti(x) = x, and T 2 (x) = 2x 2 — 1. The continuous Chebyshev 
expansion of a function u(x) on [-1,1] is 

u(x) = f) n k T k (x), (43) 

Jk=0 

where the Chebyshev coefficients, 



A 

Uk : 

= — f t i(x)T k {x)w(x)dx, 
JTCjfc J - 1 

(44) 

with the Chebyshev weight, 


w(x) = (1 - x 2 )”J, 

(45) 

and 


_ f 2 k= 0 

Cfc ( 1 otherwise 

(46) 

The orthogonality relation is 

2 

nc k 

J T k (x)Ti(x) w{x)dx = 6 kt . 

(47) 

The change of variables, 


X = CO80y 

(48) 

the definition, 


v(0) = u(co$0), 

(49) 

and (42) reduce (43) to 


OO 




v(^) = ^2 u k cosk6. 

(50) 
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Thus, the Chebyshev coefficients of u(x) are precisely the Fourier coefficients of v(0). This new 
function is automatically periodic. If u(x) is infinitely differentiable, then v{6) will be infinitely 
differentiable. Hence, straightforward integration-by-parts arguments lead to the conclusion 
that the Chebyshev coefficients of an infinitely differentiable function will decay faster than 
algebraically. Note that this holds regardless of the boundary conditions. 

Some important properties of Chebyshev polynomials are the boundary values, 



T*(±l) = (±1)*, 


(51) 


n{± i) = (±i) fc * 2 , 


(52) 

and the recursion relations, 




2x7* (z) = 

7*+i(x) + 7*_ i(i), k > 1, 


(53) 

27* (i) = 


k > 1. 

(54) 


The discrete Chebyshev expansion is based on the Gauss-Lobatto collocation points, 

*3 


x i = cos ~ N > 




The discrete Chebyshev coefficients are 


2 " 1 


kxj 


where 


N' 

_ ( 2 k = Oor N 

~\1 1 < k < N — 1 


k = 0,1,...,N, 


The inversion formula is 


N 


Uj = 53 u k coa 


k = o 


njk 
N ‘ 


The Gauss collocation points are also sometimes used: 

*u- 1 ) 


Xi — cos- 


N ’ 




(55) 

(56) 

(57) 

(58) 

(59) 


Special versions of the FFT may be used for evaluating the sums in (56) and (58). See CHQZ, 
Appendix B for a discussion of how to adapt the standard complex FFT for the Chebyshev 
transforms based on the Gauss-Lobatto and Gauss points. Schumann and Sweet (1988) discuss 
how a wide variety of sine and cosine transforms may be implemented with the standard FFT; 
this includes transforms which take advantage of special symmetries. 

The truncated Chebyshev series is defined to be 

N 

Pnu(x) = 53 UkTk(x), (60) 

i=o 

and the Chebyshev interpolating polynomial is 


N 

/jv«(i) = 53 «* ?*(*)• 


Jk=0 


(61) 
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Differentiation in Chebyshev spectral methods is based upon the relation, 

! = f>WjiW, (62) 

where the Chebyshev coefficients, tzj^ , of the derivative are related to those of the function itself 
by the recursion relation, 

2 kii k = - 4+i> k>l. (63) 

(See CHQZ, Sec. 2.4.2 for the derivation of this formula.) This relation obviously generalizes to 

2 = **-i4-i - 4+1, k > 1- (64) 

Explicit expressions for the first and second derivative coefficients are 

o oo 

4 1 = — 52 p » p > ( 65 ) 

c k 

p=k+2 
p+k odd 

and 

1 OO 

4 2) = — 52 p(p 2 - fc2 ) { *p, ( 66 ) 

c k 

p=k+2 
p+k even 


For the truncated Chebyshev series (60), the derivative coefficients, 4^> are zero f° r k — 
and the derivative coefficients for 0 < k < N — 1 may be computed in decreasing order by 


**4° = 4+2 + 2 (* + !)«*+ 1, k = N- l,N-2 0. (67) 


The Chebyshev Galerkin derivative is just d(Pptu)/dx. Unlike the Fourier case, truncation and 
differentiation do not commute: in general, d(Pffu)/dx is a polynomial of degree N — 1, whereas 
Pjq{du/dx) is a polynomial of degree N . A similar formula holds for the coefficients, 4 > f° r 
the derivatives of the Chebyshev interpolant. Higher-order derivative coefficients are computed 
by iterating on this recursion relation. 


Chebyshev collocation differentiation can also be represented by matrices, 
derivative, 


d(ijyu) 

dx 


N 


(**) = 52 D ki u i> 


1=0 


where 


f to/ffiK- 1)* + 7(**-*0 



-*»/[2(l - «?)] 

(21V 2 + l)/6 
— (21V 2 + l)/6 


0 < k, l < N, l ± k 
l<k = l<N-l 
k = l = 1 
k = l = N 


For the first 

( 68 ) 


(69) 


and for the second derivative (Ehrenstein and Peyret 1987), 


d^(+ru) 

dx 2 


{xk) = 52 D ki u i> 

1=0 


(70) 
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where 


[(-l)* +J /c*] [x\ + x k xi - 2]/[(l - xl)(x k - i/) 2 ] 1 < k < N - 1, o<l < N\ l^k 

-[(* 2 - 1)(1 - x 2 ) + 3]/[3(l - x 2 ) 2 ] l<k = l<N-l 

= | (2/3)[(-l)'/ci][(2iV 2 + 1)(1 - x«) - 6]/(l - x,) 2 k = 0, 1 < / < N 

(2/3)[(-l) w+ Vc i ][(2 J Y 2 + 1)(1 + x,) - 6]/(l + x,) 2 k = N, 0<l < N -1 
k (N 4 — 1)/ 15 k = l = 0,N 

(71) 

General expressions for the first- and second-derivative matrices for Jacobi polynomials are 
available in Deville and Mund (1989). For an explicit expression for the inverse of the Chebyshev 
first-derivative matrix (with homogeneous Dirichlet boundary conditions on one end) see Funaro 
(1988a). 

2.3 Galerkin, Tau, and Collocation Approximations 

On some domain fl with boundary dfl, consider the differential equation, 

L(u) = f, in H, 

B(u) = 0, on dCl b} (72) 

where L is a differential operator, B is a linear boundary operator, u is the solution, / is the 
data, and dQ b is that subset of dCl on which boundary conditions are enforced. The boundary 
conditions are assumed here for simplicity to be homogeneous. 

In order to use spectral methods to solve this differential problem, a set of discrete equations 
must be employed. Three strategies for deriving a discrete system are in common use: Galerkin, 
tau, and collocation. All of these may be viewed as particular cases of spectral weighted residual 
methods. The key elements are the set of trial functions, Xn, the set of test functions, IV, and 
discrete approximations, Ln and Bjy, to the operators, L and B. Another important element is 
an inner product on fi, denoted by (u,v). CHQZ, Secs. 10.3 and 10.4 have described a general 
formulation of Galerkin, tau, and collocation methods. Here we shall be content with presenting 
a straightforward prescription for these methods. 

The starting point for choosing the trial and test functions is a set, Pn, of algebraic polyno- 
mials on fi of degree less than or equal to N . Let comprise an orthogonal basis for 

and assume that these polynomials are ordered by increasing degree. For notational simplicity, 
we describe the procedure for approximations to non-periodic problems, in which case Jacobi 
polynomials are appropriate. Both the trial and test functions are contained in subsets of 
The trial functions necessarily satisfy the boundary conditions: 

Xn = {t; G Pn | Bv = 0 on 3Q&}. (73) 

The approximate solution is denoted by u N , and the residual is just 

fl(u") = L(u") - /. (74) 

The weighted residual orthogonality statement is equivalent to 

n N € X N , 

{L{u n ) -f, v) = 0, for all v G Y N . (75) 
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The Galerkin weighted residual approximation results from setting Yn = Xpj, and using the 
continuous inner product, 

(u,t;) = f uv w dCl, (76) 

J n 

where w is a weight function, in (75). Let {<j>k}k € J be a basis for X n (and in the case of 
Galerkin methods, also for Y)v). The (f>k are linear combinations of the polynomials pk such that 
the boundary conditions are satisfied. The approximate solution can be expressed as 


u N = ^4>k- 

keJ 


(77) 


Eq. (75) equivalent to the requirement that 

u N € X N , 

( L{u N Uk) = (f,4>k ), k € J. (78) 

We have assumed that Ljsf = L, as occurs for most applications. In the Galerkin approximation 
the unknowns are the expansion coefficients, u*, of the basis functions for X^. 

The tau approximation is very similar, but the test functions are not required individually 
to satisfy the boundary conditions. The set, X^ y remains as given by (73), but Yn = Pn-m> 
where M is the number of boundary conditions. Although the solution still has the form (77) 
of an expansion in terms of the trial functions, it is generally more convenient to work with it 
in the form, 


U N = Y, UkPk, 
k = 0 

B(u N ) = 0, on dn b , (79) 

of an expansion in terms of the orthogonal basis for P#. The condition (75) is equivalent to (79) 
plus 

(L(u N ),p k ) = (/,p t ), k = 0,1,. . . ,N — M. (80) 

Observe that the residual is not required to be orthogonal to the M highest-order expansion 
functions. These conditions are dropped in favor of enforcing the boundary conditions on u N . 
Note that the unknowns in the expansion are the coefficients of the basis functions for Pn and 
not of the basis functions for Xj?. In the tau method there is no need to construct explicitly the 
combinations of polynomials which satisfy the boundary conditions. 

Collocation methods are at first sight quite different. A set of collocations points Xj,j = 
1,2,..., AT, is chosen, and the approximate solution is required to satisfy (72) exactly at the 
collocation points: 


#jv(u^)(zj) = 0, at all boundary points , 

LN{u N )( x j) = f{ x j)> at internal points, (81) 

where B^(u) and Ljy(u) are the discrete approximations to B(ti) and L(u) which result from 
taking derivatives via collocation. In this approximation the unknowns are the nodal values of 
the function. 

Nevertheless, the collocation method can be placed within the framework of weighted residual 
methods. Loosely speaking, one keeps the condition (79) on the trial functions and chooses the 
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Dirac delta-functions, S(x — Xj), as the test functions. A more precise analogy is based on 
discrete inner product, 

N 

(u,v) N = J2 u i v i w J> ( 82 ) 

3=0 

which results from applying the appropriate Gaussian quadrature formula, rather than the 
continuous inner product (76). (The quadrature rule, in fact, determines the collocation points, 
Xj, and the weights, u/y.) As shown in CHQZ, Sec. 10.4.3, the collocation method can then be 
placed in the same variational framework as the Galerkin method. 

In these notes we shall customarily use u* for the expansion coefficients of Galerkin and tau 
methods (in which they are the fundamental variables), but u* for the expansion coefficients of 
collocations methods (in which they are intermediate variables used for the approximation of 
derivatives). 

2.4 Generalizations 

In many cases a mapping of the desired interval, (a, 6), onto the standard interval - (0, 2?r) for 
periodic boundary conditions and (—1,1) otherwise - is needed. A simple linear mapping is 
straightforward. However, often a nonlinear mapping can improve the quality of the approx- 
imation. Some types of mappings are discussed in CHQZ, Sec. 2.5. Of particular interest is 
the recent work of Boyd. See the references in the bibliography of CHQZ and in the present 
bibliography. 

An important unresolved issue is how severe a stretching is permitted. The current folklore 
is that the truncated series must resolve not only the function itself, but also the mapping, in 
order to achieve a good approximation. 

In more than one dimension the standard choice for the basis functions is just the tensor 
product of one-dimensional basis functions. (See, however, Bisseling and Kosloff (1988) for an 
application using a Fourier method on an isotropic grid.) This, of course, limits conventional 
spectral methods to extremely simple geometries. This limitation can be surmounted by spectral 
domain decomposition methods, which are discussed in Chapter 4. 

Some subtle issues arise in the imposition of Neumann or Robin boundary conditions and 
the treatment of coordinate singularities. See CHQZ, Ch. 3. In some circumstances there are 
advantages to enforcing a combination of the differential equation and the boundary conditions 
at the boundary. This applies to both elliptic (Funaro 1988b) and to hyperbolic (Funaro and 
Gottlieb 1988) problems. 


3 Elementary Applications to Model Problems 

In this chapter we describe the spatial discretizations that result from spectral approximations 
to simple PDE’s. For evolution problems the time-discretizations used in spectral methods are 
quite conventional; see CHQZ, Ch. 4 for a complete discussion. 


3.1 Some Fourier Methods For the Wave Equation 

An elementary example is provided by the model problem, 

du dn 


(83) 
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with periodic boundary conditions on (0,2flr), and the initial condition, 



u(x,0) = 8in(wco8z). 

(84) 

The exact solution, 

u(x,t) = 8XYl[KC08{x — i)] y 

(85) 

has the Fourier expansion, 

oo 



«0m)= £ «fc(*)«‘ fcz » 

Jb= — oo 

(86) 

where the Fourier coefficients, 

Ufc(t) = ain{^)j k {n)e~' kt , 

(87) 

and Jk{t) is the Bessel function of order k. The asymptotic properties of the Bessel functions 
imply that 

k p Uk(t) — ► 0 aa k —* oo (88) 


for all positive integers p, consistent with the result (19). As a result, the truncated Fourier 
series, 

f " 1 

«*W) = £ *k(t)e ikx , (89) 


converges faster than any finite power of 1/N. 

The prescription (78) readily extends to the semi-discrete (discrete in space, continuous in 
time) Galerkin approximation to (83) and produces the classical result, 


dxik 

— + ,kn„ = 0, 


k = 


N_ 
2 ’ 


• > 


N_ 

2 


1 . 


(90) 


Let Uj(t) denote the approximation to u(xy, t), and recall (27). Then the semi-discrete collocation 
approximation (81) to (83) is 




(xy, t) = 0, 




(91) 


N 

Fourier Galerkin 

Fourier Collocation 

Second-order 

Fourth-order 

8 

9.87 (-2) 

1.62 (-1) 

nr]o] 

9.62 (-1) 

16 

2.55 (-4) 

4.97 (-4) 

6.13 (-1) 

2.36 (-1) 

32 

1.05 (-11) 

1.03 (-11) 

1.99 (-1) 

2.67 (-2) 

64 

6.22 (-13) 1 

9.55 (-12) 

5.42 (-2) 

1.85 (-3) 


Table 1: Maximum Error for the Periodic Wave Equation 


An illustration of the superior accuracy available from the spectral method for this problem 
is provided in Table 1 and Figure 4. Shown in the table are the maximum errors at t = 1 for the 
Fourier Galerkin and the Fourier collocation methods as well as for second-order and fourth-order 
finite-difference methods. The time discretization was the classical fourth-order Runge-Kutta 
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method. In all cases the time-step was chosen so small that the temporal discretization error was 
negligible. The figure compares the second-order finite-difference and Fourier Galerkin solutions 
for N = 32 with the exact answer. Note that the major error in the finite-difference solution 
is one of phase rather than amplitude. In many problems the very low phase errors of spectral 
methods is a significant advantage. 

Because the solution is infinitely differentiable, the convergence of the spectral method on 
this problem is more rapid than any finite power of 1/N. (The error for the N = 64 spectral 
result is so small that it is swamped by the round-off error of these 60-bit calculations.) In 
most practical applications the benefit of the spectral method is not the extraordinary accuracy 
available for large N, but rather the small size of N necessary for a moderately accurate solution. 

3.2 Some Chebyshev Methods for the Wave Equation 

The PDE is again (83), but the initial condition is now 

u(x,0) = 8in(anx ), (92) 


and the boundary condition is 


u( — 1 , t) = ui,(t) = stn[ax(l - t)]. 


The function, 


u(x,t) = sin[a^r(x — t)], 

is one solution to the problem (83), (92), and (93). It has the Chebyshev expansion, 

00 

«(*»*) = X 

Jb=0 


(93) 

(94) 

(95) 


where 


Ufc(t) = — 8»«(^r - ant)J k (air). 
Ck 2 


(96) 


Clearly, the truncated series, Pjvu, converges at an exponential rate. Note that this result holds 
whether or not a is an integer. In contrast, the Fourier coefficients of u(x,f) are 


f, (t\ = * c ic*Tt 6in *(<* + k ) _ j - iart rin*(a-k) 
fcV ; 2 jt a + k 2jt a-k ' 


For non- integer a these decay extremely slowly - only as 0(k 1 ) - and the Fourier series exhibits 
a Gibbs phenomenon (see CHQZ, Sec. 2.1.4). 

For the present problem the semi-discrete Chebyshev tau formulation (79) - (80) is 

^ + 4 1) = 0. * = 0,1,. ...IV- 1, (98) 


N 

XX- 1 )* 6 * = u L{t), (99) 

Jk=0 

along with the obvious initial condition on the Chebyshev coefficients. Consistent with (80), the 
highest-order condition (k = N) for the differential equation has been dropped in (98) in favor 
of the boundary condition (99). The property (51) has been employed in (99). Note that no 
special formula is required for the derivative at j = 0 (or x = 1). In a typical finite-difference 
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method a special one-sided formula would be utilized at the right boundary. The Chebyshev 
collocation approximation is 


duj d(I N u) 


dt 


■w+ 


dx 


(xj,t) = 0, 


j = 0,1,. .. ,N - 1, 


( 100 ) 


with 


un(0 = ul(0. 


( 101 ) 


and the initial condition. 


N 

Truncated Series 

Chebyshev Collocation 

Fourier Collocation 

Second-order 

4 

1.24 (0) 

1.49 (0) 

1.85 (0) 

1.64 (0) 

8 

1.25 (-1) 

" 6.92 (-1) 

1.92 (0) 

1.73 (0) 

16 

7.03 (-6) 

1.50 (-4) 

2.27 (0) — 

1.23 (0) 

32 

1.62 (-13) 

3.45 (-11) 

2.28 (0) 

3.34 (-1) 

64 

1.79 (-13) 

9.55 (-11) 

2.27 (0) 

8.44 (-2) 


Table 2: Maximum Error for the Non-periodic Wave Equation 


Results pertaining to a = 2.5 at t = 1 for a truncated Chebyshev series, a Chebyshev col- 
location method, a Fourier collocation method, and a second-order finite-difference method are 
given in Table 2. For this non-periodic problem Fourier, spectral methods are quite inappropri- 
ate (and the failure to converge in the maximum norm reflects the Gibbs phenomenon), but the 
Chebyshev spectral method is far superior to the finite-difference method. 

The Chebyshev collocation points are the extreme points of Tn(x). Note that they are not 
evenly distributed in x, but rather are clustered near the endpoints. The smallest mesh size scales 
as 1 /N 2 , and the magnitude of the largest eigenvalue of the Chebyshev first-derivative operator 
grows as N 2 . While the Chebyshev Gauss-Lobatto distribution contributes to the quality of the 
Chebyshev approximation and permits the use of the FFT in evaluating the series, it also places 
a severe time-step limitation on explicit methods for evolution equations. 


3.3 Some Explicit Chebyshev Methods for the Heat Equation 
The heat equation, 

du __ d 2 u 
dt~ dx*' 

is the natural parabolic linear model problem. For our simple example the spatial domain is 
(-1,1), the initial condition is 

ti(x,0) = stnarx, (103) 



and the boundary conditions are 


The exact solution is then 


ti(—l, t) = u(+l,t) = 0. 


u(x,0)=e ^strurx. 

The semi-discrete tau equations are 

* = 0,l,...,J\T-2, 


du k _ (2) 
~dt~ uk ’ 


(104) 

(105) 


(106) 
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with 


(107) 


N N 

«* = X) «* = °- 

fc=0(mod 2) Jk=l(mod 2) 

(The notation A: = 0 (mod 2) means that only the contributions from the even values of k are 
included in the sum.) The Chebyshev coefficients, uj^, of the approximation to the second 
spatial derivative can be obtained from ti* by two applications of the recursion relation (63). 
In this tau approximation the dynamical equations for the two highest-order coefficients are 
dropped in favor of the equations for the boundary conditions. Eq. (107) follows from (104) 
and the property (51). The Chebyshev collocation approximation is 



d 2 (h yu) 
dx 2 


(*/.*)» 




(108) 


with 


u Q (t) = u N (t) = 0. 


(109) 


For both the tau and collocation approximations above, implementation of explicit time- 
discretization schemes is straightforward. The results at t = 1 are given in Table 3. The 
maximum errors shown there have been boosted up by the factor e* so that they represent 
relative errors. On the whole the collocation results are the best. It goes almost without saying 
that finite-difference results are far inferior to any of these spectral approximations. 


N 

Truncated Series 

Tau 

Collocation 

8 

2.44 (-4) 

1.61 (-3) 

4.58 (-4) 

10 

5.76 (-6) 

2.12 (-5) 

8.25 (-6) 

12 

9.42 (-8) 

3.19 (-7) 

1.01 (-7) 

14 

1.14 (-9) 

3.35 (-9) 

1.10 (-9) 

16 

1.05 (-11) 

8.39 (-11) 

2.09 (-11) 


Table 3: Maximum Error for Chebyshev Approximations to the Heat Equation 


3.4 An Implicit Chebyshev Method for the Heat Equation 

The time-step restriction for explicit Chebyshev methods for the heat equation is very severe, 
scaling as 1/N 4 . Implicit time-discretizations overcome this stability limit at the price of requir- 
ing the solution of implicit equations. 

We write the Chebyshev tau approximation resulting from a Crank-Nicolson discretization 
of (106) - (107) as 

- 1a* uj^(* + Ai) + u fc (* + At) = + «*(*)» k = 0,l,...,N -2, (110) 

with 

N 

x; u fc (i+A0=o, (in) 

k=0(mod 2) 

N 

u k {t+At) = 0. (112) 

Jt=l (mod 2) 
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We re-write (110) as 


(113) 


- A ii k = ft, 

where A = 2/ At, the coefficients on the left-hand side are at t + At, and 


(114) 


Eqs. (Ill) - (113) form an implicit system for the coefficients for k = 0, 1, . . . , N, with (66) 
used to express uj^ in terms of u* . The corresponding linear system is upper triangular. The 
solution process for such a system requires N 2 operations. 

A far more efficient solution procedure is obtained by re-arranging the equations. By invoking 
the recursion relation (63) twice one obtains 


Cfc-2 

4k(k - 1) 


A 

fk—2 


+ ( 1 ~ 2(J?- 1))‘‘ + 

, e *+2 / 

2(fc 2 -i) /fc + 4;b(Jb-i-i) /fc+2 ’ 


gjb+2^ 

4k(k + 1) 


«*+ 2 = 


k — 2, 3, ... , IV, 


(115) 


where 


e k = 


( 


1 

0 


0< Jfc< W-2 

k> N- 1 


(116) 


in place of (113). (See CHQZ, Sec. 5.1.2 for the complete derivation.) Note that the even and 
odd coefficients are uncoupled. By writing the equations in the order (111), and then (115) for 
k = 2,4, .. . ,N one obtains a linear system which is quasi- tridiagonal: the first row is full, and 
the remaining rows have a standard tridiagonal structure. A similar ordering is chosen for the 
odd coefficients. This ordering has been chosen to minimize the round-off errors arising from a 
specially tailored Gauss elimination procedure which performs no pivoting (and works from the 
“bottom up” rather than the more customary “top down”). Assuming that the coefficients in 
(115) have already been calculated, the solution of both systems takes only 161V operations, a 
substantial improvement over the 0(N 2 ) cost of the upper triangular system. 

The coefficient of u* in (115) is the largest coefficient and it is desirable for it to be on the 
main diagonal. The system is not diagonally dominant and, in practice, round-off errors are a 
mild problem: typically 4 digits are lost for N — 128. The accuracy may be increased through 
the use of iterative improvement or double-precision. 

The 0(N ) operation count for a single step of an implicit Chebyshev tau approximation to 
the one-dimensional heat equation is comparable to the cost of a finite-difference method. The 
transformation to physical space takes an extra 0(N log 2 N) operations. Hence, the cost of this 
spectral method is larger than that of a finite-difference method by a factor NflogiN . However, 
the spectral method requires far fewer grid-points to achieve a high accuracy. 

The class of problems for which such efficient implicit Chebyshev tau methods apply is 
quite limited. Only the simplest variable coefficients are permitted. Tuckerman (1989b) has 
provided a characterization of such problems and given a general prescription for manipulating 
the implicit equations into a banded form. 


An implicit collocation method does not reduce to a quasi-tridiagonal system. In this ap- 
proach the cost of solving the implicit equations is 0(iV 3 ). 
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3.5 Chebyshev Methods for the Poisson Equation 

As usual the discussion will begin with a linear model problem, but this time in two spatial 
dimensions. That problem is the Poisson equation, 

d 2 u d 2 u _ 
dx 2 + dy 2 ~ *' 

on the square, (— 1,1)*, with homogeneous Dirichlet boundary conditions. The choice, 

f(x,y) = — 2x 2 «inxx simry, (118) 

corresponds to the analytical solution , 

u(x,y) = sinirx siniry. (119) 



A Chebyshev expansion in both x and y is appropriate for this problem. Assume, for 
simplicity, that polynomials of degree less than or equal to N are used in both directions. The 
Chebyshev tau approximation is 


4 0) + 4 J) = /™» m,n = 0,1,..., iV — 2, 


with 


N 


N 


^2 Umn = = °) » = 0, 1, . . . , N, 

m—0(mod 2 ) m=l(mod 2 ) 

N N 

^ ^mn = ^ Umn =0, m = 0, 1, . . . , N. 

n=0{mod 2 ) n=l(mcd 2 ) 


( 120 ) 

( 121 ) 

( 122 ) 


In this notation, and denote the Chebyshev coefficients of and respectively. 

ox 2 dy 2 

The boundary conditions (121) - (122) contain 4 redundant conditions corresponding to the 4 
corner points of the square. The collocation approximation reads 

(*;> y*) + ( a i» yt) = /(»;•> Vk), j,k = l,2,...,N-l, (123) 

u jk = 0, j,k= 0 or N, (124) 

where x } - and yt denote the Chebyshev-Gauss-Lobatto points in x and y, respectively. 


N 

Truncated Series 

Tau 

Collocation 

8 

HBlIilSM 


WE5MSM 

10 

Mlflgj)— 

■■K-iajs 


12 



Msasm 

14 

NKEIIGM 


KEggEfl | 

16 

— ' ™W 

HHBH1 



Table 4: Maximum Error for Chebyshev Approximations to Poisson’s Equation 
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N 

Truncated Series 

Tau 

Collocation 

8 

2.88 (-4) 

2.79 (-3) 

117 (-4) I 

10 

6.79 (-6) 

5.26 (-5) 

2.33 (-6) 

12 

“ 1.09 (-7) 

8.86 (-7) 

3.12 (-8) 

14 

1.34 (-9) 

1.09 (-8) 

3.27 (-10) 

16 

“ 1.19 (-11) 

9.15 (-11) 

2.73 (-12) ' 


Table 4: Maximum Error for Chebyshev Approximations to Poisson’s Equation 


3.6 Fourier Methods for the Burgers Equation 

Our last, and not so elementary, example is for the Burgers equation, 

du du d 2 u 


(125) 


on (0, 2x) with periodic boundary conditions. The test problem is chosen so that the exact 
solution is 

(a ~ ct,t+ 1) 

«(*. 0 = -2 ^.?, — . < 126 ) 


where c = 1, and 


<t>(x — ct, t + 1) ’ 
<Kx,t)= f) e -[*-(2"+lW 3 /4W 


(127) 


Unlike the previous examples, this problem is nonlinear. The Fourier Galerkin approximation 
to the PDE is 


du 

Tt + 


Upiquq = -k 2 vu k , 


»— f-.f” 1 - (128) 


P+ 9 =* 


The convolution sum in this equation requires 0(jV 2 ) operations if evaluated in a straightfor- 
ward manner. Fortunately, transform methods permit it to be evaluated in only 0(N log^N) 
operations. 

Let us consider a general one-dimensional convolution sum of the form, 


E a a 

Up t;,. 

p+q=k 

_ S N i 

p>q—— — 1 


(129) 


The basic approach is to transform u k and v k to physical space, to perform there a pointwise 
multiplication, and then to transform the result back to Fourier space. We introduce the discrete 
transforms, 




Uj= £ y= 0,1,...,AT- 1, 


(130) 
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vj= £ 
*«-* 

and define the physical space product, 

1 

o' 

II 

(131) 

= U , V i> 

j = 0,l,...,N- 1, 

(132) 

and its discrete Fourier transform, 



3=0 

*=_* — — 1 
2 ’ ’ 2 

(133) 

Use of the discrete orthogonality relation (25) 

leads to 


OO 

w k = £ v, 

p+q=k 

P,9=-f f" 1 

OO 

+ £ « P V 

p+q=k±N 

(134) 


The first term on the right-hand side is the desired result. The second term is called the aliasing 
error. Hence, this naive transform method does not produce a true Galerkin method. For this 
reason it is referred to as a pseudospectral method. 

The transform method, however, can be modified to produce a true de-aliased, Galerkin 
convolution sum. Two techniques are available. One is called the 3/2-rule and it involves 
performing the transforms (130) and (131) on a grid with M > | N points. In this case one simply 
pads the original Fourier coefficients tt* and v * with zeroes for y < k < y before performing 
the transforms. The multiplication (132) at physical points is performed on M points, as is the 
transform (133). CHQZ, Sec. 3.2 demonstrate that the result of this method is a set of Fourier 
coefficients, W *, which agree with the desired coefficients, tb*, for A; = — y, ...,y — 1. This 
transform method costs roughly 50% more than the naive transform method discussed above, 
but it does produce a de-aliased result. 

The second de-aliasing technique involves applying the naive transform method twice, once 
as described above and a second time using a grid in physical space which is shifted with respect 
to the standard one by n/N. The average of these two results is alias-free. (See CHQZ, Sec. 3.2 
for the details.) 

In multiple dimensions the 3/2-rule generalizes in an obvious manner. However, the phase 
shift de-aliasing technique requires evaluations on multiple grids, e.g., on 8 separate grids in three 
dimensions. Moreover, the phase shift method is the more expensive of the two. CHQZ, Sec. 
7.2.2 furnish a detailed comparison of de-aliasing techniques for three-dimensional convolution 
sums. 


A Fourier collocation method for the Burgers equation is much more straightforward. It 
reads 


du i(*\ , ^ d 2 (I N u) 

*(v + u i aZ (*;.*) = V —Q x i (*;>*)> 


dt 


j = 0,l,...,N -1. (135) 


One can show (CHQZ, Sec. 3.2.5) that the collocation method (135) is algebraically equivalent 
to the Fourier approximation (128) employing the pseudospectral transform method. Colloca- 
tion methods do include aliasing errors. However, for problems with even more complicated 
nonlinearities than those of the Burgers equation, collocation methods are much simpler to 
implement than Galerkin methods. 
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The aliasing error can be removed from the collocation method for this problem. One simply 
forces the upper-third of the coefficients in transform space to be zero at every time step. (This 
can be achieved at nominal expense here if the viscous term is treated implicitly and a transform 
method is used for the solution.) This approach is called the 2/3-rule. 

Figure 5 shows typical results for second-order finite-difference and for Fourier collocation 
methods, and Table 5 gives a comparison between spectral and finite-difference methods. 


N 

Fourier Collocation 

Second-order 

Fourth-order 

16 

2-1 (-1) [ 

2.9 (-2) 

2-0 (-1) 

32 

~ 3.5 (-2) 

8.2 (-2) 

2.5 (-2) 

64 

3Xp) 

2.6 (-2) 

5.2 (-4) 

128 

6.1 (-8) 

3.3 (-3) 

2.5 (-5) 


Table 5: Maximum Error for Fourier Approximations to Burgers Equation 


There has been a long controversy over the seriousness of the aliasing errors. Recent numer- 
ical and theoretical results establish that aliased and de-aliased spectral methods yield the same 
asymptotic rate of convergence. CHQZ, Secs. 4.6, 7.2.4, 7.3.5, and 11.3 summarize the results 
available as of a few years ago. Recently, Maday and Quarteroni (1989) have extended these 
results to the Korteweg-de Vries equation. 


4 Domain Decomposition Methods 

The previous chapters have concentrated on spectral methods for problems in simple domains. 
There have been a number of recent developments on the use of spectral techniques in more 
general geometries and even on their use in simple geometries in ways that improve resolution 
for problems with disparate scales. The basic idea is to partition the complete domain of 
the problem into several subdomains. The approximation is spectral if increased accuracy is 
obtained by increasing the order of approximation in a fixed number of subdomains, rather 
than by resorting to a further partitioning. One situation in which this domain decomposition 
approach is useful is illustrated in Figure 6. In this particular example it is clear that at least 
two subdomains are required in order to use spectral methods at all. In the case illustrated in 
Figure 6 the partitions intersect only at the boundaries. It is also possible to use subdomains 
which overlap. This is illustrated in Figure 7. Additional advantages may arise by using separate 
subdomains instead of a single subdomain which is their union. The partitioning illustrated in 
Figure 6 leads to a distribution of Chebyshev collocation points which improves the resolution 
of the approximation. Moreover, one would expect the corresponding algebraic problem to be 
better-conditioned because there is a less extreme ratio of the largest to smallest grid spacings. 
Finally, the use of subdomains facilitates the implementation of spectral methods on parallel 
computers, especially those with local memory. 

A crucial aspect of any domain decomposition method is the manner in which solutions on 
contiguous domains are matched. Patching methods take a classical (pointwise) view of the 
differential equation. If the equation has order d, then at the interface of contiguous domains 
the solution and all its derivatives of order up to d — 1 must be continuous. For second-order 
problems this is typically enforced by requiring that the solution, ti, be continuous and that its 
normal derivatives, du/dn y match at the interface, in the sense that they are equal and opposite. 
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The condition on the normal derivative may be replaced with a continuity condition on any 
other directional (but not tangential) derivative. These continuity conditions are discretized by 
enforcing them at selected points, and thus, are satisfied exactly by any approximation. 

Alternatively, the differential equation may be posed variationally. The standard weak for- 
mulation is based on an integration- by- parts procedure over the entire physical domain. This 
implicitly assumes continuity of u and its first d — 1 derivatives at all points, and, in particular, 
at the interfaces between subdomains. In most cases continuity of u is built into the choice 
of trial functions. Continuity of the derivatives of u occurs as a natural interface condition. 
By reversing the integration-by-parts procedure on the decomposed domain, one can deduce 
directly from the variational principle the continuity of the derivatives of u at the interfaces. In 
a variational domain decomposition method, the natural interface conditions are satisfied with 
increasing precision as the order of approximation is increased. 

In an overlapping domain method only the continuity of the function is imposed. But this 
matching occurs on the boundary of the intersection between two adjacent domains. In the 
discrete case this is enforced at selected points on the boundary. 

4.1 Patching Methods 

A one-dimensional domain decomposition for the interval, H = (a, 6), illustrated in Figure 8. It 
is broken into two subdomains, Qi = (ai,a 2 ) and CI 2 — (< 12 , as), where ai = a and as = 6. The 
solution to a given differential equation on H is denoted by u. Its restriction to Cl 8 is denoted 
by u a , for s = 1,2. The integer, iV,, is used to indicate the degree of approximation on Q 8 . The 
approximate solution on □ is denoted by u N and its restriction to Cl 8 by u?. 

Consider a patching method for the linear problem, 


d?u du 

Lu = -t'-r-T + Of— + 
dx l dx 

Au = /, in H. 

(136) 

The ODE problem may be formulated as 



Lui = f, 

ui(oi) = 0, 

in fix, 

(137) 

Lu 2 = f, 
u 2 (o 3 ) = 0, 

in fl 2 , 

(138) 

Ui(o 2 ) = 
dux , 

-^M) = 

u 2 (a 2 ), 
du 2f s 

(139) 


One can show that this formulation is equivalent to (136). The Chebyshev collocation equations 
for (137) - (139) are apparent. They constitute N\ + N 2 + 2 independent conditions for the 
Ni + N 2 + 2 unknown nodal values. 

Let us consider now a nonlinear generalization, namely, 

= /, (140) 

where the flux, G(u), is given by 

G(u) = g(u) - (141) 


| 


24 


and </(u) is a nonlinear function of ti. The appropriate interface conditions are continuity of u 
and G. Thus, the condition, 

= G(u 2 (a 2 )), (142) 

replaces the last condition in (139). 

An alternative patching method consists of replacing the pointwise flux balance condition 
(142) by an integral version over the two subdomains adjacent to the interface point: 

G(u i(oi))-/ /(ui)dx = G(u 2 (a 3 )) + [ f{u 2 )dx. (143) 

J a\ J aa 

The correct discrete enforcement of (143) requires that appropriate quadrature formulas be 
employed. Gottlieb and Streett (in preparation) have studied this matter carefully. 

Additional considerations arise in two or more dimensions. The interface is now a curve 
rather than a point. The collocation points on the two domains may not coincide, as illustrated 
in Figure 9. In this case appropriate interpolation formulas are needed. Moreover, there may 
be internal cross points, as illustrated in Figure 10. Special care is needed at these points. The 
integral flux balance condition (143) extends easily to such points. Further details are provided 
in CHQZ, Sec. 13.2. 

Hyperbolic problems require different interface conditions than elliptic ones. Consider the 
one-dimensional wave equation (83) on fi. A patching method for this problem reads 


with 


dui dui 

_+_ = 0 , onn u 

tti(<*i,i) = u L {t), 

(144) 

du 2 du 2 

HT + 'ax =0, <mn! ’ 

(145) 

ui(a 2 ,t) = u 2 (o 2 ,t). 

(146) 


The collocation equations are obvious, except at the domain interface. One choice is to average 
the contributions from fii and Q 2 - However, the temporal stability is dependent upon the 
particular values of N\ and A 2 . The preferred choice is to collocate (144) rather than (145) 
since information is travelling to the right; in other words, the interface condition is upwinded. 


The upwind interface condition extends in a natural way to quasilinear hyperbolic systems, 
such as 

£ + A(u)£ = 0, on fi, (147) 

where u E R n , and A(u) is an n x n matrix, by a diagonalization procedure. Convergence 
of the multi-domain discretization of general hyperbolic systems of conservation laws, including 
the case in which the interface coincides with a shock, is discussed in Quarteroni (1989). An 
iteration-by-subdomain procedure, taking into account the direction of transfer of characteristic 
information at the interface, is also discussed there. 

Patching methods for the spatial part of hyperbolic problems have generally been coupled 
with explicit time-discretizations. Elliptic problem, of course, require implicit methods as do 
implicit time-discretizations of parabolic problems. The matrices which represent the global 
algebraic system have a block structure due to the domain decomposition, and only adj acent 
subdomains (or blocks) are coupled. Block Gauss elimination, then, should be considered. So, 
too, should static condensation (Przmieniecki 1963). These methods are most attractive in the 
context of time-dependent problems in which the factorization of the global matrices need only 
be performed once. 
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An influence-matrix technique is a useful alternative. This enables the solution to the global 
system to be obtained at the price of two local solutions to Dirichlet problems for each sub- 
domain. The influence-matrix is constructed by solving on each subdomain Dirichlet problems 
that have boundary data which is 1 at a particular interface point and 0 elsewhere, and then 
computing the resultant errors in the patching condition at all the interface points. In the first 
step, arbitrary values are assigned to the interior interface nodes. Then the residuals for the 
patching conditions are evaluated. These will be nonzero, but the influence- matrix is used to 
infer the correct interface values. The local solutions are then re-computed with the correct 
interface values and the desired solution is obtained. 

CHQZ, Sec. 13.2.3 cover some of the early solution algorithms. More recent work on the 
influence-matrix method has been reported by Pulicani (1986), Gottlieb and Hirsh (1988), and 
Phillips and Karageorghis (1989). Block preconditioners of the influence-matrix have been 
discussed by Quarteroni and Sacchi-Landriani (1988, 1989). 

As mentioned before, domain decomposition methods are often used to improve resolution 
in problems possessing both large- and small-scale structures. An example of this (taken from 
Macaraeg and Streett 1986) is the solution of the Burgers equation, written in the form, 


du ld(ti 2 ) 
dt + 2 dx 

d 2 u 

x € (-1,1), 

(148) 

«(-M) = 

u(M) = o, 

t > 0, 

(149) 

tx(x,0) = — stn7rx, 

x G (-1,1), 

(150) 


and with v = 0.01/x. Table 6 presents a comparison of a conventional single-domain Chebyshev 
collocation method with a computation which employed three domains. (A mapping was also 
employed to improve resolution of the internal “shock” .) Observe that the multi-domain solution 
requires only half as many points to achieve accuracy comparable to the single-domain method. 
Moreover, the decay rate of the multi-domain solution exhibits spectral accuracy. 


Method 

Discretization 

Relative Error in Maximum Slope 

single-domain 

64 

1.31 (-4) 

multi-domain 

12/13/12 

" 1.99 (-4) 

multi-domain 

20/21/20 

3.23 (-5) 

multi-domain 

32/33/32 

2.14 (-6) 


Table 6: Comparison of Single-domain and Multi-domain Solutions to the Burgers Equation 


4.2 Variational Methods 

The essential aspects of variational spectral domain decomposition methods can be gleaned from 
an examination of their application to the one-dimensional Helmholtz problem, 




(152) 


An equivalent, variational formulation of this problem is that u is the solution to 

[ . dudv + \ uv \£ x — f f v dx, forallveE , 

J n dx dx Jq 

where E is the space of all functions which vanish on 3Q and which, together with their first 
derivatives, are square integrable over O. On each subdomain the approximation is a polynomial 
of degree less than or equal to N a . 

In a variational domain decomposition method the space of trial (and test) functions Xpf is 
given by 

X N = {v E C°( U)\v 8 G for s = l,2 and v = 0 on dCl} f (153) 

where v s denotes the restriction of v to fi*. The space, X^ 7 is therefore composed of continuous 
functions which are piecewise polynomials defined on the decomposition of the physical domain, 
fl. Eq. (152) is simply replaced by its discrete counterpart in which the test and trial spaces 
are Xp? rather than the full space E: 

+ \u N v)dx = (f,v), forallv€X N , (154) 

where (/, v ) is a convenient approximation to / n fvdx, usually obtained by numerical quadrature. 

This is a variational formulation with trial and test functions which are continuous across 
element (or subdomain) boundaries. Flux continuity, i.e., continuity of the normal derivative, 
at element interfaces, is not satisfied for fixed N, but only as part of the convergence process, 
i.e., when N tends to infinity. 

In the spectral-element method, the discrete approximation satisfies (154), and, on each sub- 
domain (here called element), u N has a polynomial expansion of high degree. In the spectral- 
element method, the basis in is formed by Lagrange interpolants using the Legendre collo- 
cation points. 

After a global node numbering and an assembly of the elemental equations, the algebraic 
representation of (154) takes the form, 


CU = BF. (155) 

The unknown vector, U , contains the values of the discrete solution at all grid-points except the 
end points of the interval. The matrices, C and B, are positive-definite, symmetric and banded, 
the bandwidth being determined by the largest N,. 

When a two-dimensional problem is considered, the Cartesian products of the Legendre- 
Lobatto nodes are used on each subdomain. Then Lagrangian interpolants at these points can 
be represented using tensor products of Legendre polynomials. Spectral-element methods have 
also been developed which use isoparametric elements. More details are provided in CHQZ, Sec. 
13.3.2. 

CHQZ, Sec. 13.5.2 demonstrate an equivalence between certain patching methods and vari- 
ational methods. The equivalent patching condition in variational methods follows from appro- 
priate quadrature rules. It now appears clear that patching methods yield the best results when 
the interface condition is derived variationally. Funaro (1988b) presents some recent results on 
this subject. 

The articles by Maday and Patera (1989) and Fischer, Ho, Karniadakis, Ronquist, and Patera 
(1988) provide recent summaries of the latest developments. Among these is the construction 
of nonconforming spectral-element methods by Maday, Mavriplis, and Patera (1988). 
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Early implementations of the spectral-element method favored direct methods, in particular, 
the static condensation method, for the solution of the implicit system represented by (155). 
More recently, iterative methods have been employed. See the discussion in Sec. 5.2.5 and the 
article by Fischer et al (1988). 

4.3 Alternating Schwarz Method 

This class of domain decomposition methods uses overlapped domains, as illustrated in Figure 
8. For elliptic problems an iterative procedure is clearly required in order for the solution on 
the individual subdomains to converge to the restriction of the solution on the entire domain. 
CHQZ, Sec. 13.4 discuss some of the techniques currently in use. 

Recently, Kopriva (1989a) presented a spectral domain decomposition method on a combi- 
nation of patched and overlapped grids for hyperbolic equations. Figure 11 illustrates the grid 
that he employed for a two-dimensional hyperbolic system, and Figure 12 demonstrates the 
achievement of spectral accuracy in the steady-state. Kopriva (1989a) discusses the merits of 
two different procedures for handling the overlapped grids. 

4.4 Recent Developments 

Spectral domain decomposition methods have been an extremely active area of research in the 
last few years. A recent collection of papers has appeared in a special issue of Appl. Niimer. 
Math, (cited individually in the bibliography). Here we summarize some of the recent work 
which was not covered in CHQZ, Ch. 13. 

Several methods have been developed in which finite-difference or finite-element methods are 
used in some domains, whereas spectral methods are used in others. Kopriva (1989b) coupled 
upwind finite-difference methods to improve resolution in the vicinity of a shock wave with 
spectral methods in the smooth part of the flow. Bernardi, Debit, and Maday (1987) proposed 
and analyzed a method which used finite-elements near a singularity and spectral methods 
elsewhere. 

Gastaldi and Quarteroni (1989) devised a well-posed mathematical formulation of a problem 
in which the equation itself differed on the subdomains, in particular hyperbolic in one domain 
and parabolic in the other. They also analyzed a spectral domain decomposition method for 
this application. 
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5 Solution of Implicit Equations 


The solution of implicit equations is an important component of many spectral algorithms. For 
steady problems this task is unavoidable, while spectral algorithms for many unsteady problems 
are only feasible if they incorporate implicit (or semi-implicit) time-discretizations. We concen- 
trate on linear systems, assuming that nonlinear ones are attacked by standard linearization 
techniques. 

We focus on the elliptic equation, 

Au — A u = /, (156) 

where / is a function of x and A is a constant. The Poisson equation corresponds to A = 0, and 

an implicit time-discretization of the heat equation gives rise to such a system with A > 0 (see 

Sec. 3.4 for the one-dimensional case). 

If Fourier series are used for any spatial direction, then a Fourier transform in that direction 
eliminates the derivative and simply alters A. For this reason, we discuss only Chebyshev ap- 
proximations here. However, CHQZ, Sec. 5.1 discuss direct methods for Fourier approximations 
at some length, including cases involving separable, variable coefficients and mappings. 

Spectral collocation approximations lead to a linear system, 

LU = F, (157) 

where U and F are vectors consisting of the grid-point values of u, / and any boundary data, 
and L is a matrix (constructed as a tensor prpdp^t matrix in two or more spatial dimensions). 
A similar linear system is obtained for Galerkin and tau approximations, but now U and F 
are vectors consisting of the expansion coefficients of u, / and the boundary data, and L is the 
appropriate matrix in transform space. 

The spectral linear systems are essentially full. The matrices for one-dimensional problems 
are dense and the matrices for multi-dimensional problems are tensor products of these full 
matrices. As Deville and Mund (1989) have stressed, there are a large number of zeros in the 
spectral matrices for multi-dimensional problems; however, the bandwidth is nearly maximal. 
Gaussian elimination may, in principle, be applied, but it requires 0(JV 3d ) operations (where 
addition, subtraction, multiplication, and division are counted as separate and equal operations) 
and 0(N d ) storage, where d is the dimension of the problem. (We assume, for simplicity, that 
the number of degrees of freedom in each spatial dimension is N.) In the first section of this 
chapter we discuss some direct solution techniques which can yield the solution to (157) in 
0(N d ) y 0(N d log2N) y or at worst 0(N d+1 ) operations with at most 0(N d ) additional storage. 
The remainder of the chapter is devoted to a discussion of iterative techniques. These require 
0(N d log2N ) operations per iteration and 0(N d ) additional storage. 

5.1 Direct Methods 

Our objectives in this section are to explain the principles underlying the basic direct solution 
techniques. We shall call a solution “efficient” if it enables the solution to (157) to be obtained 
in at most 0(N d log2N) operations. This makes the cost of solving (157) comparable, even for N 
large, to the cost of typical explicit spectral operations such as differentiation and the evaluation 
of convolution sums. In many cases, a solution cost of 0(N d + l ) is still acceptable in the sense 
that it only overwhelms the cost of other spectral operations for values of N of 128 or so. 

An important consideration is whether only a few or else a large number of solutions to 
(157) with different data, F y are sought. The latter case is typical of semi-implicit methods 
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for unsteady problems: hundreds, or even thousands of solutions to a single linear system are 
required. In such situations, it is reasonable to invest a substantial amount of calculations on a 
pre-processing stage that greatly reduces the subsequent cost of solving (157). 

Let us consider (156) on a square with homogeneous Dirichlet boundary conditions. The 
collocation approximation to this can be written 

AU + UB - XU = F, (158) 

where U is the (N - 1) x (N - 1) matrix (tiy) for t, j = 1, 2, . . . , N- 1, F is defined similarly, A is 
the second-derivative operator (in x) in which the boundary conditions have been incorporated, 
and B is the transpose of the second-derivative operator (in y). 

Systems of the form (158) are solvable by Schur-decomposition. An orthogonal transforma- 
tion is used to reduce A to block-lower-triangular form with blocks of size at most 2. Similarly, 
B is reduced to block-upper-triangular form. If P and Q denote the respective orthogonal 
transformations, then (158) is equivalent to 

A'U' + U'B' - XU' = F' t (159) 

where 

A' = P t AP, 

B' = Q t BQ, 

U' = P t UQ, 

F 1 = P t FQ. (160) 

The solution process has 4 steps: 1) reduction of A and B to real Schur form (and de- 
termination of P and Q); 2) construction of F' via (160); 3) solution of (159) for U'\ and 4) 
transformation of U 1 to U via the inverse of (160). 

The first step can be accomplished via the QR algorithm in approximately 40 N 3 opera- 
tions. Step 3 requires 2 N 3 operations, and steps 2 and 4 take 4 N* operations apiece. A single 
solution requires roughly 50 N s operations. Hence, step 1 is the most time-consuming. When 
the same problem must be solved repeatedly, then step (1) need only be performed once, in a 
pre-processing stage. The matrices A',B',P, and Q may then be stored and used as needed. In 
this case a complete solution takes 10N S operations. 

The matrix-diagonalization approach is similar to the Schur-decomposition method. The 
difference is that the matrices, A and B, in (158) are diagonalized rather than merely reduced 
to block-triangular form. An algebraic problem of the form (159) is obtained with (160) replaced 
by 


A! = P~ l AP = \ A) 

B 1 = Q~ l BQ = A B) 

U' = P~ X UQ, 

F' = P~ l FQ , (161) 

where Ax is the diagonal matrix with the eigenvalues of A on the diagonal. Thus, we have 

A A U' + U' K b - XU' = F 1 . (162) 

The matrices, P and Q, are not necessarily orthogonal and their columns consist of the eigen- 
vectors of A and B, respectively. 
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The matrix-diagonalization scheme for (158) consists of the same 4 steps as the Schur- 
decomposition method, except that the first, pre-processing stage also requires that the eigen- 
vectors and the inverse transformations be computed. This takes an additional 8 N 3 operations. 
Step 3 takes only 3JV 2 operations since the system is diagonal, and steps 2 and 4 require 4 N 3 
operations apiece, as before. 

For collocation problems requiring multiple solutions, the matrix-diagonalization method has 
the advantage of taking only 80% of the solution time of the Schur-decomposition method: 8 N* 
operations. Moreover, the entire solution process - steps 2, 3 and 4 - is extremely simple and 
can be optimized readily. In fact, most computer libraries contain assembly-language routines 
for all the requisite operations. The third stage of the Schur-decomposition method is more 
complicated. 

In the case of tau approximations to (156), further gains in efficiency are possible. The 
discrete problem may be written in the form (158) where U is the (N — 1) X ( N — 1) matrix, 
(u nm ), consisting of the Chebyshev coefficients of u (minus those used to enforce the boundary 
conditions). The matrix, F, is defined similarly, and A and B are the representations in trans- 
form space of the second-derivative operator (with the boundary conditions used to eliminate 
the 2 highest-order coefficients in each direction). 

In the case of Dirichlet (or Neumann) boundary conditions, the even and odd modes decouple. 
Thus, A,J5,P, and Q contain alternating zero and non-zero elements. This property may be 
exploited to reduce the cost of both the pre-processing step (by a factor of 4) and the matrix 
multiplies (by a factor of 2). The cost of steps 2 through 4 is thus 4 AT 3 . 

The cost of the solution stages may be halved again by performing the diagonalization in 
only 1 direction and resorting to a standard tau solution in the other. Thus, (158) is reduced to 

AU’ + U'G - XU 1 = F\ (163) 

where 

U' = UQ , 

F* = FQ y (164) 

instead of (161). The system (163) decouples into N — 1 systems of the form (113). Each of 
these may be reduced to a system like (115) and solved accordingly in 16AT operations. The 
cost of the solution process is essentially halved, to 27V' 3 operations, since the number of matrix 
multiplies is cut in two. 

In these algorithms, as indeed with matrix computations in general, the accumulation of 
round-off error is a concern. Haidvogel and Zang (1979) reported the loss of 3 to 4 digits (for 
N between 16 and 64) with the Schur-decomposition method. These were recovered through 
iterative improvement. Since the computation of eigenvectors can be a sensitive process, double- 
precision is advisable for the pre-processing stage of the matrix-diagonalization method. 

Both methods can be generalized. The use of Neumann or Robin boundary conditions is 
straightforward. However, with Robin boundary conditions the even and odd modes do not 
de-couple; hence, some of the economies of the tau method are lost. These methods can also be 
applied to separable equations of the form, 

£(*>£)+£(*<)-*•=/• 

As it happens, these methods are more attractive in three-dimensional problems than in two- 
dimensional ones. Suppose that the number of degrees of freedom in each direction is N . The 
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pre-processing cost is some large multiple of JV 8 . In two dimensions, the solution cost is a small 
multiple of N s > and typical explicit spectral calculations take 0 (N 2 log 2 N) operations. Thus, the 
pre-processing cost is substantially larger than the cost of a single solution. In three dimensions, 
the pre-processing cost remains 0(A 3 ), the solution cost is a small multiple of N 4 and typical 
explicit spectral calculations require 0(N z log2N) operations. Thus, the pre-processing cost may 
even be smaller than the cost of the solution phase. Similarly, the extra memory required for P 
and its inverse is proportionately smaller in three dimensions than in two. 

5.2 Iterative Methods 

The direct methods discussion above have limited application compared to iterative methods. 
Invariably, iterative solution schemes are applied to spectral methods employing collocation. 


5.2.1 Richardson Iteration 


The fundamentals of iterative methods are perhaps easiest to grasp for the simple model problem, 


d * 2 


(166) 


on (0, 27r) with periodic boundary conditions. The Fourier collocation approximation to the 
left-hand side of (166) at the collocation points is 


r— 1 


E pV’’''- 


.= -? + ■ 


The spectral approximation to (166) may be expressed as 

L U = F, 


(167) 


(168) 


where U — (uo, uj, . . . , ti^v_i), F = (/o> /i, • • • , /n~i), and L represents the Fourier spectral 
approximation to —d 2 /dx 2 . 

A Richardson’s iterative scheme for solving (168) is 

yn+l = yn + _ L yn ) ? (jgg) 

where V n is the n th iterate and a; is a relaxation parameter. The eigenfunctions of L are 

Mp) = ^" ,k , p=-f,..-,f -1, (170) 

with the corresponding eigenvalues, 


A P -P 2 , 


N N , 

p ~ 2 ’ 2 lf 


(171) 


The index p has a natural interpretation as the frequency of the eigenfunction. 

The error at any stage of the iterative process is V n - U; it can be resolved into an expansion 
in terms of the eigenvectors of L. The error obeys the relation, 

(^"+1 - u) = G[V n - U ), (172) 
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where the iteration matrix of the Richardson scheme is given by 


G — I - uiL. 


(173) 


The iterative scheme is convergent if the spectral radius p of G is less than 1. 

Each iteration reduces the p th error component to i/(A p ) times its previous value, where 

i/(A) = 1 - u;A. (174) 


(This is referred to as the damping factor.) The optimal choice of u results from minimizing 
I u(X) I for A e [A mtn , A ma z], where A min = 1 and A mai = ^ . (One need not worry about 
the p = 0 eigenfunction since it corresponds to the mean level of the solution, which is at one’s 
disposal for this problem.) The optimal relaxation parameter for this procedure is 


oj = 


: + A r 


It produces the spectral radius, 


P = 


: + A f 


K — 1 
K + 1 ’ 


(175) 


(176) 


where 


k — 


Amaz 

Amtn 


(177) 


is referred to as the spectral condition number. Unfortunately, k ~ and p - 1 “ W’ which 
implies that 0[N 2 ) iterations are required to achieve convergence. 


5.2.2 Preconditioning 

The primary cause of the inefficiency of the straightforward Richardson method is the rapid 
increase with N of the spectral condition number. This can be alleviated by “preconditioning” 
the problem, in effect solving 

{H-'MViU = ( H~ 1 M)F (178) 

rather than (168). A preconditioned version of (169) is 

V n+ i =V n + u>H~ 1 M(F - LV n ). (179) 

In practice the inverse of the preconditioning matrix, H , is never explicitly required; instead one 
solves 

H(V n+1 - V n ) = u>M(F - LV n ), (180) 

The effective iteration matrix is 

G = I - uH-'ML. (181) 

One obvious requirement for H is that (180) can be solved inexpensively, i.e., in fewer 
operations than are required to evaluate LV ; it is also desirable for MV to be inexpensive to 
evaluate. The second requirement on the preconditioning matrices is that H~ 1 M be a good 
approximation to L -1 , i.e., that the spectral condition number of H^ML be small. 

A finite-difference approximation to L is one useful choice for the preconditioner. In this case 
H represents a finite-difference approximation to the left-hand side of (166) and M is the identity 
matrix. A finite-element discretization is another, and usually superior, choice. In practical 
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terms the principal difference between finite-difference and finite-element preconditioning is that 
the latter includes a weighting of the residuals. Here H is the appropriate stiffness matrix and 
M is the mass matrix. For one-dimensional problems, only 0(N ) operations are required to 
solve (180) and MV takes only O(N) operations. The evaluation of LV takes 0(N log 2 N) 
operations. Hence, the preconditioning work is a minor part of the cost of a single iteration. 

Specific formulas for finite-difference preconditioners are fairly straightforward and are dis- 
cussed in CHQZ, Sec. 5.2. Canuto and Pietra (1987) and Deville and Mund (1989) both provide 
thorough discussions of the appropriate formulation of the finite-element preconditioners. 

Simple estimates can be made of the effectiveness of various preconditioners for the model 
problem (178) since the eigenfunctions of the preconditioned problem have the form e tkx >\ Let 
an eigenvalue of the preconditioned operator, H~ 1 ML ) be denoted by A. For example, a simple 
calculation indicates that the eigenvalues arising from second-order finite-differences are 

Ap = (pAi/2) 2 /s»n 2 (pAi/2), p = -y, ■ ■ . , y - 1. (182) 

Clearly, A mfn = 1, and A maz = tt 2 /4. Table 7 summarizes the properties of second-order finite- 
difference (FD2), fourth-order finite-difference (FD4), and linear finite-element (FE1) precondi- 
tioning. It turns out that these simple estimates are also good predictors of the performance of 
these preconditioning techniques for Chebyshev approximations to non-periodic problems. Note 
that the finite-element preconditioning achieves the smallest condition number, just 1.44. 


Preconditioner 

Amm 

A max 

P 

FD2 

1.00 

2.47 

0.424 

FD4 

1.00 

1.85 

0.298 

FE1 

.693 

1.00 

0.181 


Table 7: Properties of Preconditionings for the Model Problem 


Ronquist and Patera (1987b), in the context of the spectral element method, have used 
scaled Jacobi preconditioning, i.e., H = a D, where D is the matrix containing the diagonal 
elements of L, and a is a scaling parameter chosen so that the spectral radius of H matches 
that of L. The primary application of this has been in the context of multigrid methods and it 
is discussed further in Section 5.2.5. 

Central-difference preconditioners work well for problems which are dominated by even-order 
derivatives. When significant odd-order derivatives are present, however, other preconditioners 
must be employed. The Fourier first-derivative operator, preconditioned by second-order central 
differences, has the eigenvalues 

N N 

A P = (pAx)/sin{pAx), p l. (183) 

Now, A m , n = 1, but h-max — oo. No iterative scheme can succeed for this spectrum. Two 
alternatives are to employ a staggered grid, or to apply the preconditioning on a grid with more 
points, typically 50% more, than the basic grid. The former alternative is discussed in CHQZ, 
Sec. 5.2.2. 

In the latter case, effective preconditioning of second derivatives in a mixed-order problem 
is retained. Figure 13, taken from Streett and Macaraeg (1986), shows how the eigenvalue 
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spectrum of the preconditioned Chebyshev first-derivative operator is compressed in shifting 
the central-difference preconditioning grid from the original one to ones which are 20% and 50% 
finer. 


5.2.3 Conventional Iterative Methods 

Much better iterative methods than the elementary Richardson scheme discussed above are 
available. Both non-stationary Richardson and Chebyshev acceleration of Richardson iteration 
are superior: in effect the term k in (176) is replaced by y/n. Variational iterative schemes such 
as minimum residual and conjugate gradient are other possibilities. They have the advantage of 
being parameter-free. See CHQZ, Sec. 5.3 for a thorough discussion. 

5.2.4 Multi-dimensional Problems 
Consider now the generalized Helmholtz problem, 

V • (fl(x)Vti) - A u = /, (184) 

in two-dimensions. The discrete approximation to (184) may also be written in the form (168), 
with obvious notation. Assuming that N collocation points are employed in each direction, the 
evaluation of LV takes 0(N 2 log2N) operations, and the term MV costs 0(N 2 ) operations. The 
direct solution of (180) takes 0(A 3 ) operations by a banded matrix solver after the 0(N 4 ) cost 
of factoring the matrix has been paid. In this case the cost of a single iteration is dominated 
by the preconditioning. The situation is even worse in three dimensions, in which the cost of 
evaluating LV is only 0{N z log2 N) y the cost of solving (180) is 0(N 4 ), while the factorization 
expense climbs to 0(N 5 ). 

Of course, the finite-difference or finite-element preconditioner may be inverted approxi- 
mately by an iterative scheme, preferably conjugate gradient or multigrid. The essential ap- 
proach is the iterative solution of the preconditioned system (180). Experience indicates that 
the iterative solution of the preconditioned system need not be carried out to a high degree of 
convergence; typically a reduction in the residual of the preconditioned problem by two to three 
orders of magnitude suffices for convergence of the overall iteration. Multigrid solution of the 
finite-difference preconditioned system is particularly effective. 

Alternatively, in two and three dimensions the economics may favor trading off less effective 
iterative schemes (accepting larger spectral radii) in favor of less expensive preconditioners. 
Applications to date in this area resort to starting with a finite-difference approximation and 
using incomplete-LU decompositions, ADI relaxation, or alternating line relaxation (see CHQZ, 
Sec. 5.4). 

ADI and alternating line relaxation are well-suited to vector and shared-memory parallel 
processors; incomplete-LU decompositions are not. Of course, these methods would surely ben- 
efit from replacing the finite-difference approximation with a finite-element one. An interesting 
alternative, which has apparently not yet been tried for spectral methods, is to employ the poly- 
nomial preconditioning strategy of Dubois, Greenbaum, and Rodrique (1979). Their approach 
is to write H = D(I — G), where D is diagonal, and to use a polynomial in G as an approxi- 
mation to (/ — G) -1 . In this case the trivial inversion of D and the approximate inversion of 
(J--G) both vectorize and parallelize well. This approach does require ||G|| < 1 for convergence. 
Wong (1989) showed that for solving fourth-order discretizations, an effective strategy is to use 
the second-order finite-difference approximation to determine G. This use of the second-order 
finite-difference approximation will not succeed as a preconditioner for spectral discretizations, 
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since the condition number of H X ML is greater than 2. (It must be less than 2 to satisfy the 
requirement on G.) However, finite-element approximations should work well (see Table 7). 


5.2.5 Multigrid Methods 


Iterative schemes for spectral collocation equations can be accelerated dramatically by apply- 
ing multigrid concepts. Briefly put, multigrid methods take advantage of a property shared by 
a wide variety of relaxation schemes - potential efficient reduction of the high-frequency error 
components, but unavoidable slow reduction of the low-frequency components. This slow con- 
vergence is the outcome of balancing the damping of the lowest-frequency eigenfunction with 
that of the highest-frequency one in the minimax problem described after (174). The multigrid 
approach takes advantage of the fact that the low-frequency modes (|p| < N / 4 for a coarse grid 
with only half as many points as the fine grid) can be represented just as well on coarser grids. It 
settles for balancing the middle- frequency eigenfunction (|p| = N j 4) with the highest-frequency 
one (jp| = N/ 2), and hence damps effectively only those modes which cannot be resolved on 
coarser grids. In (175) and (176), A min is replaced with A mt( j = A(^). The optimal relaxation 
parameter in this context is 

U MG = T • (185) 

A max i A mid 

The multigrid smoothing factor, 




Amax ^ 


mid 


A max “f“ A mi *d 


(186) 


measures the damping rate of the high-frequency modes. In this example, hmg = 0.60, in- 
dependent of N . The price of this effective damping of the high-frequency errors is that the 
low-frequency errors are hardly damped at all. Table 8 compares the single-grid and multigrid 
damping factors (see (174)) for N = 64. However, on a grid with N/2 collocation points, the 
modes for | p | € [^, are now the high-frequency ones. They get damped on this grid. Still 
coarser grids can be used until relaxations are so cheap that one can afford to damp all the 
remaining modes, or even to solve the discrete equations exactly. For the case illustrated in 
Table 8, the high-frequency error reduction in the multigrid context is roughly 250 times as fast 
as the single-grid reduction for N = 64. 

As in most multigrid methods the smoothing rates are distinctly slower as the dimensionality 
of the problem increases. For this simple model problem, the smoothing rate is 0.33 in one 
dimension, 0.60 in two dimensions, and 0.78 in three dimensions. 

Let us consider just the interplay between two grids. A general, nonlinear fine-grid problem 
can be written 

L*{U 7) = F f . (187) 

The shift to the coarse grid occurs after the fine-grid approximation, V* , has been sufficiently 
smoothed by the relaxation process, i.e., after the high-frequency content of the error, W — U? > 
has been sufficiently reduced. The related coarse-grid problem is 


L c {U e ) = F c , (188) 

where 

F c = R[F f - L f (V')] + L e (RV f ). (189) 

The restriction operator, i?, interpolates a function from the fine grid to the coarse grid. The 
coars6-grid operator and solution are denoted by L e and £7 C , respectively. After an adequate 
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p 

Single-Grid 

Multigrid 

1 

.9980 

.9984 

2 

.9922 

.9938 

4 

.9688 

.9750 

8 

.8751 

.9000 

12 

.7190 

.7750 

16 

.5005 

.6000 

20 

.2195 

.3750 

24 

.1239 

.1000 

28 

.5298 

.2250 

32 

.9980 

.6000 


Table 8: Damping Factors for N = 64 


approximation, V c , to the coarse-grid problem has been obtained, the fine-grid approximation 
is corrected via 

V f <- V s + P{V C - RV f ). (190) 

The prolongation operator, P, interpolates a function from the coarse grid to the fine grid. A 
common strategy for shifting between grids is illustrated in Figure 14. This is called a V-cycle. 
Usually a coarse grid has only half as many points (in each coordinate direction) as the fine grid. 

A complete multigrid algorithm requires specific choices of the interpolation operators, the 
coarse-grid operators, the relaxation schemes, and a strategy for shifting between grids. Gener- 
ally, the effective spectral radius, which measures the average reduction in residual per fine-grid 
relaxation, is somewhat slower than the smoothing factor, hmg j due to errors introduced in the 
interpolation stages. These issues are discussed at length in CHQZ, Sec. 5.5 for both Fourier 
and Chebyshev multigrid methods. 

The isotropic, fully periodic Poisson problem (A = 0) can be solved extremely efficiently by 
spectral multigrid methods. The discretization of problem (184) is called isotropic if 


Ax 2 Ay 2 Az 2 ' 

Using a residual weighting technique, the effective spectral radius can be as low as 0.11 in 
two dimensions and 0.18 in three dimensions. (Note the slower convergence with increased 
dimensionality.) This is achieved with purely explicit relaxation. Preconditioning is entirely 
unnecessary. Details are given in Brandt, Fulton, and Taylor (1985), Erlebacher, Zang, and 
Hussaini (1987), and Heinrichs (1988b). 

It is important to recognize that these effective spectral radii presume that the coefficient 
a(x) is reasonably well-behaved. As Erlebacher, Zang, and Hussaini (1987) demonstrated, if 
this coefficient is highly oscillatory, then the convergence deteriorates, and if the coefficient is 
basically random, the spectral multigrid method fails. Of course, similar behavior occurs for 
finite-difference multigrid methods for such problems. Nevertheless, numerical examples have 
been presented of successful spectral multigrid solutions of the isotropic problem with fairly 
strongly variable coefficients. 

If the fully periodic problem is anisotropic, or if the coefficient A is non-zero, then a different 
strategy must be employed. The best procedure demonstrated to date in two dimensions is to 
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employ a finite-difference preconditioner together with alternating line relaxation. Effective con- 
vergence rates between 0.3 and 0.4 have been reported by Brandt, Fulton, and Taylor (1985) and 
Heinrichs (1988b). Note that the cost of this type of preconditioner is insignificant compared 
with the cost of the spectral residual evaluation. As a preconditioner for a single-grid itera- 
tive scheme, alternating line relaxation performs poorly, since the eigenvalues associated with 
the smooth components of the solution tend rapidly to zero with N . Combining alternative 
line relaxation with the spectral multigrid scheme provides both an effective (and nearly grid- 
independent) convergence rate with an iterative scheme costing only O^N^log^N) operations 
per iteration. 

For problems with at least one non-periodic direction, preconditioning is again essential. 
Heinrichs (1988a) has demonstrated that in two dimensions alternating line relaxation produces 
the same effective spectral radius (between 0.3 and 0.4) that it achieves for the fully periodic, 
anisotropic problem. He also observed that for Dirichlet boundary conditions, minimum residual 
relaxation was more robust than Richardson relaxation. For Neumann boundary conditions, 
however, nonstationary Richardson must be used for minimum residual fails to converge. 

Spectral multigrid methods for the Legendre-based spectral-element method have recently 
been developed (Ronquist 1988, Ronquist and Patera 1987b). To date they have been based on 
the Jacobi preconditioning discussed in Section 5.2.2. In one dimension the effective spectral 
radius for the simple Poisson problem is 0.7, independent of both if, the number of elements, and 
of N , the degree of approximation in each element. In two dimensions, the multigrid condition 
number is again independent of if but it grows as N. (In contrast, this apparently grows no faster 
than N 1 / 8 for the spectral multigrid methods utilizing more sophisticated preconditioning.) The 
best two-dimensional spectral radii that have yet been achieved are above 0.9, and these have 
been for N no larger than 12. To date, this method has not been tested for variable coefficients 
anywhere near as extreme as those for which conventional spectral multigrid methods have 
succeeded. 

The essential problem is that as N grows the aspect ratio of the grid formed by the collo- 
cation points increases. This is akin to having an anisotropic problem on a uniform grid. For 
finite-difference and finite-element multigrid, as well as for the conventional spectral multigrid 
methods as discussed above, the only known effective relaxation procedures for two-dimensional 
anisotropic problems involve incomplete-LU decomposition or line relaxation. It is our opin- 
ion that the proper route to acceptable convergence rates for two-dimensional spectral-element 
multigrid involves replacing the diagonal preconditioning used to date with more effective pre- 
conditioners. 

Maday and Munoz (1989) have furnished a rigorous proof of the convergence of this spectral 
multigrid method. However, it does presume sufficient smoothness for the coefficient a, so that 
the potential failure of the method to converge for highly oscillatory coefficients is not addressed. 
(Note that even if the coefficient is well-resolved on the fine grid, it may be too poorly resolved 
on a much coarser grid for convergence.) 

As yet, spectral multigrid methods have not been applied to three-dimensional anisotropic 
problems. However, past experience with finite-difference approximations to such problems 
suggests that in three dimensions alternating line relaxation alone will not suffice, regardless of 
whether one is dealing with spectral or spectral-element multigrid; instead, alternating plane 
relaxation is needed. The article by Dendy (1987) contains a summary of some recent experience 
with finite-difference multigrid methods for anisotropic three-dimensional problems. The planar 
relaxation itself can be accomplished by alternating line relaxation within a multigrid context 
or by preconditioned conjugate gradient iteration. 

We also recommend that finite-difference preconditioning for spectral multigrid methods be 
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replaced by finite-element preconditioning. Table 7 makes it clear that superior convergence 
rates will be achieved. In fact, we project that ultimately the combination of alternating line 
(or, in three-dimensions, plane) relaxation with finite-element preconditioning will produce ef- 
fective spectral radii sis good as 0.25, even for anisotropic three-dimensional problems, but with 
reasonably well-behaved coefficients. 

6 Euler Applications 

Spectral methods have met with mixed success for problems in inviscid flow. One form of the 
Euler equations is 


f t +V-(p u) = 0, 

^M + V.(uu) + V P =0, 

+ u • Vp + 7 pV • u = 0, (192) 

where p is the density, u the velocity, p the pressure, and 7 the ratio of specific heats (an ideal 
gas law is presumed here). The overwhelming difficulty is that this is a nonlinear hyperbolic 
system, and it admits non-smooth solutions - shock waves and contact discontinuities. A priori, 
one expects spectral methods to perform poorly for such problems. Two strategies for dealing 
with shock-wave discontinuities have been employed - shock-capturing and shock-fitting. 

6.1 Shock-capturing 

In this approach the discontinuity is located in the interior of the domain. Finite-difference 
methods require some form of explicit or implicit artificial viscosity to handle the shock. Shock- 
capturing spectral methods have resorted to a combination of explicit artificial viscosity and 
sophisticated filtering techniques. 

Filtering methods for spectral methods have been investigated for linear model problems, 
and CHQZ, Secs. 8.3 and 12.1.4 summarize this work. Some recent contributions to artificial 
viscosities for spectral approximations to the periodic conservation laws have been made by 
Tadmor (1989) and by Maday and Tadmor (1989). Cai, Gottlieb, and Shu (1989) have employed 
discontinuity subtraction techniques to construct spectral schemes with bounded total variation. 
In all three cases exponential convergence for smooth solutions has been proven. Of course, 
the real challenge is to produce exponential convergence in smooth regions for discontinuous 
solutions. Thus far, numerical tests have been disappointing. 

It is as true today as it was 4 years ago (Hussaini, Kopriva, Salas, and Zang (1985a)) that 
although spectral shock-capturing methods have succeeded in producing solutions to the Euler 
equations, they have yet to demonstrate that they are superior or even comparable in efficiency 
to modern upwind finite-difference, finite-element, and finite-volume approaches. 

For the simpler inviscid problem of compressible potential flow, however, the dependent 
variable itself is continuous. In this case spectral shock-capturing methods (combined with a 
spectral multigrid solution algorithm) have proven most efficacious. See CHQZ, Sec. 8.5.1. 

An interesting alternative is the use of a domain decomposition method in which the region 
surrounding the shock is treated by an upwind finite-difference scheme and the remaining part of 
the flow by a Chebyshev spectral method. Kopriva (1989b) has demonstrated that this approach 
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produces accuracy in the spectral domains that is superior to the accuracy achievable by a pure 
upwind scheme. 

6.2 Shock-fitting 

In the spectral shock-fitting approach the burden of resolving the discontinuity is removed by 
making the shock one of the boundaries of the domain. The numerical method, then, computes 
the shape and motion of the shock front as part of the solution process. This technique is 
described in CHQZ, Sec. 8.6. 

Recently, Kopriva, Hussaini, and Zang (in preparation) in a sequel to Hussaini, Kopriva, 
Salas, and Zang (1985b), have improved the shock-fitting method used in their earlier work. 
Improvements have been made in both the shock-fitting itself and in the boundary conditions. 
They have re-visited the blunt-body problem that was reported in Hussaini, Kopriva, Salas, and 
Zang (1985b). The boundary conditions on the body were changed to a compatibility condition, 
and the complete (rather than approximate) compatibility condition was applied at the shock. 
As a result of these improvements, not only is a true steady-state achievable (to within machine 
precision), but also spectral accuracy is achieved, even without periodic filtering of the solution. 
Figure 15 displays the decay of the error in the stagnation pressure (which is known exactly) as 
the grid is resolved for Mach 4 flow past a circular cylinder. See this forthcoming work for the 
details. 

7 Navier-Stokes Applications 


The convection form of the incompressible Navier-Stokes equations is 

(193) 

(194) 
u 7 and 

Stokes 

(195) 

(196) 

The nonlinear term is no longer explicitly present. Spectral discretizations of the Navier-Stokes 
equations usually include an explicit treatment of the nonlinear term, and one may consider 
that it has been absorbed into f . For most purposes the basic algorithms can be stated for the 
Stokes problem. 

The numerical analysis of spectral methods for the steady Stokes problem, 

— uAn + Vp = f, (197) 

V • u = 0, (198) 

has been the first step in an analysis of the unsteady Navier-Stokes equations. The second step 
is the analysis of the steady Navier-Stokes equations. The inclusion of the nonlinear terms in 
the analysis is usually accomplished by resorting to an implicit function theorem (see CHQZ, 
Sec. 11.3). Analysis of fully-discrete approximations to the unsteady problems is much more 
difficult, and results are only now beginning to appear (Sacchi-Landriani and Vandeven 1989). 


ciu 

— + u • Vu — v Au + Vp = f, 
ot 

V • u = 0. 

The velocity, u = (tii,u y ,u z ), the static pressure is denoted by p, kinematic viscosity by 
a possible source term is denoted by f . 

A simpler problem, and one which exhibits most of the numerical difficulties is the 
problem, 

du 

— - i>Au + Vp = f, 

V • u = 0. 


i 


40 


7.1 Fourier Methods for Fully Periodic Problems 

Homogeneous, isotropic turbulence is perhaps the one fluid dynamical problem for which strictly 
periodic boundary conditions in all spatial directions are justifiable. Hence, Fourier spectral 
methods are ideally suited for this class of problems. Moreover, since the nonlinearities of the 
Navier Stokes equations are at worst quadratic, Fourier Galerkin methods are the most natural 
and efficient spectral techniques for this problem. 

7.1.1 Galerkin Methods 

The primitive variables are expanded in a three-dimensional Fourier series: 


u(x,t) = H"(*) e ’ kx > 

k 

(199) 

p( x >*) = Z)M0 e ’ kx » 

(200) 


k 


where k = (k x , k v , k z ), and each component ranges between — y and y — 1. The semi-discrete 
Galerkin approximation to (193) - (194) is 

+ j/fc 2 u k + tkpjj = r k , (201) 

tic • u k = 0, (202) 

where fc 2 = |k| 2 , and 

r k =-(urvu) k + f k . (203) 

The incompressibility condition implies that the pressure term in (201) is 

= ( 204 ) 

The only complication in this Orszag-Patterson algorithm is the evaluation of the nonlinear 
term, (u* Vu)^. This is done most efficiently by transform methods and costs 0(N 3 log2N) 
operations. The one-dimensional version of transform methods was discussed in Section 3.6. If 
the naive (aliased) transform method is applied, then the algorithm is a pseudospectral one. If 
the de-aliased transform method is used, then the algorithm is a true spectral Galerkin method. 

The Galerkin method can be implemented with as few as 8 three-dimensional FFT’s per 
step. A linear coordinate transformation has been developed that permits simulation of flows 
with constant strain, shear, and rotation within the confines of periodic boundary conditions. 
See CHQZ, Sec. 7.2 for more details. 

7.1.2 Collocation Methods and the Skew-symmetric Form 

Fourier collocation approximations to this problem are also possible. For these approximations 
use of the appropriate form of the Navier-Stokes equations is crucial. 

The standard convection form is given by (193), and the divergence form differs from it 
by the replacement of the nonlinear term by V • (u u). The rotation form of the momentum 
equation, given by 

^ + wxu— i/Au + V(i|u| 2 +p) =f, (205) 
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where the vorticity, 


cj = V x u, (206) 

has long been the preferred form. However, recently Zang (1989), following the work of Horiuti 
(1987), has demonstrated that the skew-symmetric form, given by 

^7 + ^u- Vu+ iv- (uu) -uAu + Vp = f, (207) 

ot l i 

is superior. The term skew-symmetric is used because the operator |v * Vu + • (vu) (for 

fixed v satisfying V • v = 0) is skew-symmetric. 

As discussed in Chapter 2, for Fourier Galerkin methods the projection (spatial discretiza- 
tion) and differentiation operators commute, but for collocation methods they do not. This 
means that for the semi-discrete momentum equation, the Fourier Galerkin convection, diver- 
gence, rotation, and skew-symmetric forms are algebraically identical. On the other hand, for 
collocation methods none of these forms are algebraically identical. 

Consider now the conservation properties for the ideal inviscid {y = 0), fully periodic 
case. The ideal Navier-Stokes equations conserve linear momentum, f udV, and kinetic en- 
ergy, / 5 -|u| 2 dV\ F° r the semi-discrete ideal equations, a Galerkin method conserves the discrete 
counterparts of both quantities. However, for a collocation method the semi-discrete conser- 
vation properties depend upon the exact form of the momentum equation which is discretized. 
In particular, for the rotation and skew-symmetric forms, both quantities are conserved, for 
the divergence form only linear momentum is conserved, and for the convection form neither 
quantity is conserved. 

The conservation of kinetic energy is virtually mandatory for a simulation to be numerically 
stable in time. Thus, for quite practical reasons use of the convection or divergence forms is 
ruled out. Although either of these forms alone is numerically unstable, a method which uses 
the convection and divergence forms on alternate time steps appears to be well-behaved. We 
shall refer to this version as the alternating form. 

Only 6 derivatives are required for the evaluation of the nonlinear terms for a collocation 
approximation in the rotation form, whereas 18 derivatives are needed for the skew-symmetric 
form. The convection, divergence, and alternating forms take 9 derivatives. (By invoking the 
incompressibility constraint, the number of derivatives for the convection and skew-symmetric 
forms can be reduced by 1.) Thus, the rotation form takes the least work for the evaluation of the 
nonlinear terms, the alternating form slightly more, and the skew-symmetric form appreciably 
more. Recall, however, that there is additional work required for the pressure and viscous terms. 

Zang (1989) illustrated the difference between the Fourier collocation methods using the 
rotation and skew-symmetric (actually the alternating) forms, and compared them with the 
Galerkin method on a homogeneous turbulence problem involving a mean flow with uniform 
shear. The initial condition consisted of a random divergence-free velocity field with a prescribed 
three-dimensional energy spectrum. The difference between the algorithms was dramatically 
evident in the temporal evolution of the turbulence intensity, 

, = <u| +«;+«;>*, (208) 

where (•) denotes a spatial average. Figure 16 displays the results for the rotation, alternating, 
and Galerkin forms of the nonlinear term for the test problem on 64 s , 96 s , and 128 s grids. The 
rotation form on the 64 s grid proved incapable of sustaining the turbulence at the proper level. 
The alternating form also suffered a slight loss of turbulence intensity, but a much milder loss 
than the rotation form. The same trend is apparent on the 96 s grid, but the actual differences are 
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less. The rotation form is still unsatisfactory, but now the alternating form is quite acceptable. 
On the 128 3 grid, even the rotation form now appears adequate. Clearly, all the (stable) forms 
of the momentum equation converge to the same result. However, the absolute error level is 
much larger for the (collocation) rotation form. 

In the presence of physical boundaries, the ideal conservation laws still hold with the bound- 
ary fluxes taken into account. However, for the requisite expansions in Chebyshev polynomials 
the Galerkin projection and differentiation operators do not commute, and a semi-discrete con- 
servation law, even accounting for the boundary fluxes, does not hold precisely (see CHQZ, Sec. 
4.5). Nevertheless, it holds to a sufficient degree of precision to warrant a preference for the 
conservation form of the semi-discrete equations. Zang (1989) demonstrated that the advantages 
of the skew-symmetric form extend to methods employing Chebyshev polynomials. Ronquist 
(1988) showed that it has advantages as well in spectral-element methods. In both these appli- 
cations the CPU time is dominated by the solution of the implicit equations, so that the use of 
the skew-symmetric form in place of the rotation form adds a negligible cost to the calculation. 


7.2 Fourier-Chebyshev Methods for Partially Periodic Problems 


No-slip boundary conditions are required at walls. In directions normal to the walls, then, the 
flow is non-periodic and Chebyshev expansions must be used. Customarily, the diffusion term 
in non-periodic directions is treated implicitly to overcome an otherwise prohibitive stability 
restriction on the time-step. In this section we consider problems with walls in only one direction 
and periodic boundary conditions in the other directions. Some important applications are to 
plane and circular Poiseuille flow, to Taylor-Couette flow with periodic conditions in the axial 
direction, and to the parallel boundary layer. 

The primary difficulty of algorithms for wall-bounded incompressible flows is the simultane- 
ous enforcement of the incompressibility constraint and the no-slip boundary condition. Three 
approaches are in common use: splitting methods, coupled methods, and Petrov-Galerkin meth- 
ods. 

The incompressibility constraint is most easily but least rigorously satisfied in splitting meth- 
ods. Let the spatial domain be denoted by Cl and the walls by 30. Modern splitting methods 
take the form of an advection-diffusion step, 


C7 U 

— — h u • Vu — l/Au 
ot 

= f, in Q, 

(209) 

U = g, 

on dfi, 

(210) 

followed by a pressure correction step, 



du 

+ Vj> = 0, 

in H, 

(211) 

V • u = 0, 

in Cl, 

(212) 

u • n = 0, 

on dCl, 

(213) 

(|? + Vp).f = 0 , 

on dCl, 

(214) 

where n and t denote unit vectors normal to and tangential to the wall, respectively. (No-slip 
boundary conditions on the tangential velocity may not be imposed in the pressure correction 
step; doing so produces an ill-posed problem.) At the end of the full step the velocity will be 
divergence-free, but there will be a slip velocity at the walls. 
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If the intermediate boundary conditions are simply g = 0, then this slip velocity will be 
O(At). The boundary conditions (210) on the intermediate step can be chosen to reduce the 
final slip velocity to 0(At 2 ), via 

g • n = 0, 

g • f = A t Vp • f, (215) 

or to 0(At 5 ) via 

g * n = 0, 

g • f = (At Vp + At 2 |-Vp) • f . (216) 

at 

(These formulas presume that a backward Euler time-discretization is applied to (211) and 

(214).) 

The big advantage of these splitting techniques is that they require the solution of only Pois- 
son equations (for the pressure) or Helmholtz equations (from a Crank-Nicolson discretization of 
the viscous term). These positive-definite, scalar equations are much easier to solve numerically 
than the indefinite, coupled equations that arise in unsplit methods. The implicit equations are 
diagonal in the periodic directions after a Fourier transform is applied, so that if there is only 
a single non-periodic direction, then one must solve just a set of independent one-dimensional 
problems. 

If no-slip boundary conditions are imposed on the intermediate step, i.e. g = 0, instead 
of using (215) or (216), then the O(At) slip velocity corresponds to splitting errors which are 
O(l) near the boundary for the normal pressure gradient and diffusion terms. They appear 
to cause no serious errors in the channel flow problem. However, the boundary errors produce 
serious inaccuracies in Taylor-Couette flow, which is shear, rather than pressure-gradient, driven. 
No such difficulties have yet surfaced when appropriate intermediate boundary conditions are 
employed. See CHQZ, Sec. 7.3.2 for further discussion. 

One way to avoid the splitting errors is to integrate the incompressible Navier-Stokes equa- 
tions in a single step that couples the divergence-free constraint with the momentum equations. 
The numerical difficulty of coupled methods is that one must invert a larger set of equations (it 
involves the pressure as well as the three velocity components), which is indefinite. The influence- 
matrix technique enables this solution to be achieved in the constant viscosity case by a series 
of scalar Poisson and Helmholtz solutions. This is discussed in great detail for channel flow in 
CHQZ, Secs. 7.3.1 and 11.3.3, and a recent general formulation has been provided by Tuckerman 
(1989a). Ehrenstein (1988) has derived the eigenvalue spectrum of the influence-matrix method 
for the two-dimensional streamfunction-vorticity formulation of the Stokes equations, and he 
has used these results to analyze the stability of the fully discrete approximation. Lacroix, 
Peyret, and Pulicani (1987) have used the influence-matrix method in conjunction with domain 
decomposition in streamfunction-vorticity variables. The Uzawa scheme, which is discussed in 
Sec. 7.3, is another option for avoiding splitting errors. 

Many of the numerical problems caused by the incompressibility constraint can be avoided by 
an expansion in functions which are divergence-free. Such schemes are characterized as Petrov- 
Galerkin methods since both the trial and test functions satisfy the boundary conditions, but 
they are different sets of functions. Thus far, spectral Petrov-Galerkin methods have been 
applied to pipe flow, plane channel flow, curved channel flow, and the parallel boundary layer. 
The efficiency of these methods is dependent upon the bandwidth of the matrices which arise 
from the implicit treatment of the viscous terms. In the examples cited above, the bandwidth 
is quite small, roughly of order 10. This requirement has dictated the use of special Jacobi 
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polynomials rather than Chebyshev ones in pipe and boundary layer flow. As a consequence, 
transform methods are not applicable in the non-periodic direction. Moreover, in even slightly 
more general cases, the matrices can be completely full. A description of this method is available 
in CHQZ, Sec. 7.3.3, and its rigorous numerical analysis (for the steady problem) is presented 
by Pasquarelli, Quarteroni, and Sacchi-Landriani (1987). 

7.3 Chebyshev Methods for Fully Non-periodic Problems 

Spectral methods for problems which are non-periodic in more than one direction tend to be 
more expensive than the methods discussed previously in this chapter due to the added expense 
of solving implicit equations with multiple non-periodic directions. (Recall the discussion in 
Chapter 5.) For simplicity the discussion in this section will focus on two-dimensional problems. 

The first consideration in formulating such a method is the choice of the collocation mesh. 
The various options are illustrated in Figure 17. A non-staggered mesh, in which both velocity 
components as well as the pressure are defined at the Gauss-Lobatto points, given by (55), is of 
course the most straightforward to implement. However, such a discretization admits spurious 
pressure modes, and generally requires some form of pressure boundary condition. Spurious 
modes for the pressure are characterized as those non-constant pressures which have vanishing 
x-derivatives at all the interior nodes for u and vanishing y-derivatives at all the interior nodes 
for v. Such modes have no effect upon the velocity since the interior velocity nodes are where 
the momentum equations are enforced. 

Schemes with velocities collocated on the usual Gauss-Lobatto points, but with the pressure 
defined on the Gauss points, given by (59), are said to use a half-staggered mesh. Although this 
discretization circumvents the requirement for pressure boundary conditions, it does admit one 
spurious mode, and has seen little use. 

In a fully-staggered mesh, the pressure is collocated on the Gauss points, whereas the velocity 
components use a mesh which is a tensor product of the usual Gauss-Lobatto points and the 
extended Gauss points. The latter consist of the Gauss points (59) plus the boundary points 
-1 and +1. The u velocity component is defined for the Gauss-Lobatto points in x and the 
extended Gauss points in y, and the v component is defined for the reverse combination. CHQZ, 
Sec. 7.4.1 discuss one method for carrying out differentiation on the extended Gauss points. 
Although cumbersome to implement, the fully-staggered mesh has the advantage of eliminating 
spurious modes. 

Spurious modes have been of more concern to theorists than to practitioners of numerical 
methods for incompressible flow. They are distasteful to theorists for they prevent the usual 
type of error bounds to be placed on the approximate pressure. Nevertheless, the admission 
of spurious modes need not disqualify a numerical scheme. All that’s required is that they be 
filtered from the pressure. Morchoisne and his colleagues have identified the spurious modes for 
both two- and three-dimensional problems and have described how they may be filtered from 
the solution. Bernardi, Canuto, Maday, and Metivet (1988) provide a recent discussion of both 
the practical details of the filtering procedure and the rigorous numerical analysis of spectral 
methods on non-staggered grids. Additional discussion is available in CHQZ, Secs. 7.4.1 and 
11.3.1. 

Splitting methods, particularly on a non-staggered mesh, are straightforward extensions of 
the algorithm discussed in Sec. 7.2. The equations (211) -(214) for the pressure correction step 
may be combined to form a Poisson equation for the pressure, 

Ap=i-V.u ) in U, (217) 
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where u is understood to be the velocity resulting from the first, advection-diffusion step. (Al- 
though we have written this Poisson equation in continuous form, it should actually be con- 
structed directly from the discrete forms of (211) - (214). CHQZ, Sec. 7.3.2 describe this 
process in detail in the context of a partially periodic problem.) Once the pressure is found 
from (217), the velocity is updated via (211), (213), and (214). For the non-staggered mesh, 
(217) requires a pressure boundary condition. Streett and Hussaini (1989) have used the ho- 
mogeneous Neumann condition. This is justified on the grounds of consistency with the normal 
component of (211) at the boundary. 

Coupled methods, of course, avoid splitting errors. The influence-matrix technique is one 
approach to solving the resulting implicit equations. See Tuckerman (1989a) for a thorough dis- 
cussion of influence-matrix methods for both tau and collocation discretizations of the primitive 
variable formulation, and see Ehrenstein and Peyret (1987) for such methods in streamfunction- 
vorticity variables. Influence-matrix methods are subject to an ambiguity concerning the treat- 
ment of the points at the corners of the domain. This issue is discussed in the above two 
references. 

An alternative coupled approach has been applied to spectral and spectral domain decompo- 
sition methods recently and is particularly efficient for the steady problem. The discrete version 
of (197) - (198) may be written 

-CQ + 9P = F, (218) 

OQ = 0, (219) 

where C denotes v times the discrete vector Laplacian, Q the discrete gradient operator, D the 
discrete divergence operator, and Q = (U,V). Solving (218) for Q and substituting into (219) 
yields 

(D£- 1 g)P= {DZ~')F. ( 220 ) 

This general approach is known as the Uzawa strategy. 

The solution scheme for (220) is motivated by the observation that V and Q are first-order 
operators, whereas £, is second-order. The Uzawa operator, DL ~ x £, thus is essentially zeroth- 
order, aside from the effect of boundary conditions which are enforced in the operator £. As a 
result the eigenvalues of the Uzawa operator are tightly clustered, and the condition number of 
the problem is quite small. Apparently, the eigenvalues of the discrete Uzawa operator in (220) 
are clustered more tightly for spectral discretizations than for the finite-element discretizations 
for which the algorithm was originally devised. 

Conjugate gradient relaxation is appropriate for the iterative solution of (220); it yields rapid 
convergence for problems with limited eigenvalue spread. (Actually, since the operator is not 
symmetric the relaxation scheme is more appropriately referred to as a truncated conjugate gra- 
dient iteration.) Of course, for each of these conjugate gradient iterations, two scalar Helmholtz 
problems must be inverted. The typical spectral radius observed for the conjugate gradient 
iteration in spectral applications of the steady Uzawa algorithm is between .05 and .15, with 
only an extremely weak degradation with mesh refinement. 

Streett and Hussaini (1989) have employed the spectral Uzawa algorithm on a fully stag- 
gered mesh to solve for steady states in three-dimensional nonlinear Taylor-Couette flow. This 
problem has one periodic and two non-periodic directions. The nonlinear convection term is 
incorporated by a Newton iteration, with each step involving the solution of a Stokes problem. 
By applying a Fourier transform in the periodic direction, the Stokes problem de-couples into a 
set of independent two-dimensional problems. See Streett and Hussaini (1989) for the details of 
the use of the fully staggered mesh. The Uzawa algorithm has also been employed in numerous 
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spectral-element calculations of the steady, incompressible Navier-Stokes equations (Maday and 
Patera 1989). 

Unfortunately, for the unsteady Stokes problem (195) - (196) the Uzawa algorithm is much 
less efficient. Eq. (218) becomes 

-CC-A I)Q + gP = F y ( 221 ) 

where A ~ 1/At, and (220) generalizes to 

[p{c - xiy 1 g\p = [p(c - xiy^F. ( 222 ) 

For moderate to small time steps (At <?C 1), the unsteady Uzawa operator, P(C - A/) _l £, is 
essentially second-order. Consequently, its eigenvalues are widely spread; in fact, its condition 
number scales as N 4 . In practice, the truncated conjugate gradient method fails to converge for 
all but the smallest problems. 

Streett and Hussaini (1989) developed a multigrid technique to alleviate this difficulty. Eq. 
(221) is the basic fine-grid equation, which in multigrid notation is written as (187). Their 
relaxation scheme was conjugate gradient iteration, and standard multigrid methodology was 
used to take advantage of rapid smoothing on coarse grids - see Section 5.2.5. Two relaxations 
were typically performed at each level of the fixed V-cycle; relaxations were performed only 
during the grid-coarsening phase of the cycle. (In terms of Figure 14, Nj = 2 and N u — 0.) 
Natural spectral interpolation operators were used for restriction and prolongation. On the 
coarsest mesh, a direct solution of the Uzawa operator was used; the operator was generated, 
inverted and stored in a preprocessing step. 

The multigrid/conjugate gradient algorithm is quite robust. For solution of the steady 
problem, the multigrid algorithm requires about the same amount of machine time as the single- 
grid/conjugate gradient algorithm; the improvement in convergence rate just offsets the overhead 
of the algorithm. For the unsteady problem, a small degradation occurs in the convergence rate 
of the multigrid method with decreasing A t relative to the steady case, but the machine time 
to reach a fixed residual level is roughly constant over a wide range of time-steps. (Recall that 
the single-grid scheme is not a viable option at all for the unsteady problem.) This is because 
the initial residual levels are lower for the case of smaller At (with the solution at the previous 
time-step used for the initial guess). 

The staggered-mesh multigrid/conjugate gradient coupled scheme is considerably more ex- 
pensive per time-step than the non-staggered splitting algorithm. On equivalent grids, the 
coupled scheme takes about 50-100 times longer per time-step than the split scheme. The split 
scheme is much more efficient per time-step because it requires only the solution of 11 scalar 
uncoupled Helmholtz/ Poisson equations, and these can be solved very rapidly by the matrix- 
diagonalization method (see Sec. 5.1). However, the coupled scheme is somewhat more accurate 
than the split scheme, due to the splitting errors and imposed pressure boundary conditions of 
the latter. To reduce the magnitude of these errors in the split scheme, the time-step chosen 
must be between 10 and 100 times smaller than that required purely by accuracy of the tem- 
poral discretization. Thus neither scheme is clearly superior; comparison of these schemes for 
transition simulations is ongoing. 

One spectral domain decomposition application of particular interest incorporates a non- 
reflecting outflow boundary treatment (Streett and Macaraeg 1989). The lack of viable non- 
reflecting outflow boundary conditions for the Navier-Stokes equations, ones capable of benignly 
passing large unsteady disturbances out of the computational domain, i.e., without disrupting 
the flow upstream of the outflow boundary, has limited the progress that large-scale numerical 
simulations have made in investigations of transition in shear flows. While the knowledge gained 
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through temporal simulations has been significant (see Zang and Hussaini 1987 for a recent re- 
view), the streamwise periodicity assumptions underlying these simulations limit their usefulness 
as precise, quantitative, predictive tools for spatially-developing flows. Moreover, the neglect of 
non-parallel effects in temporal simulations prevents them from addressing such fundamentally 
spatial effects as receptivity and roughness. 

The outflow boundary treatment of Streett and Macaraeg (1989) first breaks the computa- 
tional domain into two regions, the “primary” domain, which is the region of physical interest, 
and a “buffer” domain which is appended onto the downstream boundary of the primary do- 
main. This scheme is motivated by the recognition that, for incompressible flow, the spatial 
ellipticity of the Navier-Stokes equations has two sources - the viscous terms and the pressure 
field. The upstream influence influence of streamwise diffusion is ameliorated by multiplying 
the streamwise viscous term in the buffer domain by a function which smoothly reduces this 
term to zero by the end of this domain. The viscous term remains unmodified in the primary 
domain. Additionally, the convective velocity in the nonlinear advection term, which in the 
primary domain is the sum of mean and disturbance velocities, is smoothly modified in the 
buffer domain so that only the (strictly out-going) mean velocity appears at the buffer domain 
outflow boundary. The upstream influence of the pressure field is suppressed by attenuating 
the forcing function of the pressure Poisson equation in the buffer domain. A homogeneous 
Neumann boundary condition is utilized at outflow. The interface between the primary domain 
and the buffer domain is treated in a standard domain decomposition manner. 

Although a matrix-diagonalization direct solver can be used to invert the discrete Helmholtz 
equations resulting from the usual viscous terms, inversion of the modified viscous operator in 
the buffer domain apparently must be performed iteratively. (The eigenvalues of the modified 
viscous operator are complex and some of them are nearly zero.) Preconditioned minimum 
residual iteration is used, with the appropriate second-order central finite-difference operator 
used as a preconditioner. This operator is inverted using LSOR. Both iteration schemes converge 
rapidly, due in part to the scalar term of 0(1/ At) from the discrete time derivative yielding 
overwhelming diagonal dominance of the system. 

Disturbances from inconsistencies in pressure boundary conditions are generally confined to 
narrow regions near the boundary; in this application, the region of disturbance falls well inside 
the buffer domain. 

This outflow boundary condition has been verified against the predictions of spatial, linear 
stability theory for plane Poiseuille flow in a channel with the streamwise length of the primary 
domain 20 times the channel half-height; the buffer domain is the same length. The Reynolds 
number for the simulation is 1500 based on channel half-height. The numerical simulation em- 
ployed 33 points in the streamwise direction for the primary domain and 21 points in the normal 
direction. At the inflow boundary a periodic disturbance is imposed which is the least stable 
eigenmode of the corresponding linear stability problem; the maximum velocity perturbation of 
this inflow condition is roughly 0.01% of the mean center-fine velocity. The temporal period 
of the disturbance is about 19.25. Figures 18 and 19 illustrate the performance of the outflow 
technique. Figure 18 displays contours of disturbance vorticity at equally-spaced time intervals. 
The leading wave of the disturbance clearly transmits from the primary domain into the buffer 
domain without noticeable upstream disruption. Additionally, the leading wave has reached the 
end of the buffer domain by the last frame of the sequence. Figure 19 shows a plot of surface 
disturbance vorticity, taken at about 24 cycles after the start of the simulation. The discretiza- 
tion used for the primary domain in this computation was 45 X 33. The results of the simulation 
are compared in the figure against the “exact” result from linear theory. As can be seen, the 
velocity of propagation and the amplitude decay are predicted well, except perhaps near the 
endpoints of the plot, or the corners of the domain. 
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7.4 Compressible Problems 

Spectral numerical simulations of compressible, viscous flow have begun to emerge. Erlebacher 
and Hussaini (1989) studied transition on a flat plate boundary layer in supersonic flow. They 
used a straightforward, fully explicit scheme. Gauthier (1988) and Dang and Loisel (1987) 
have developed semi-implicit spectral collocation algorithms for compressible flow and have 
applied them to two- and three-dimensional problems in compressible convection. Simulation of 
homogeneous compressible turbulence (CHQZ. Sec. 8.4) are becoming quite common. 

Another method currently being examined for the solution of the time-dependent compress- 
ible Navier-Stokes equations is the fully-implicit technique of Streett and Macaraeg (in prepa- 
ration). The main idea here is to formulate the time advancement in a backward-Euler form, 
then to solve the resulting equations iteratively. The equation set we wish to solve at any time 
level is ^ 

A U n+1 + L(U n+l ) = 53 PkU n ~ k , (223) 

k = 0 

where the /?*’ s are coefficients of the (p+ l) tfc -order backward time difference. Currently, second- 
to fifth-order differencing is used. Eq. (223) is solved by first Newton-linearizing about a 
provisional value, U* : 

[0iU + A{U')\V m - X) PkU n ~ k = R m , (224) 

k—0 

then driving R m to zero using preconditioned iteration. The provisional value U* is updated 
several times to reduce linearization errors. When this process converges, the value V m is 
taken as U n+1 . Since the operator, A, is dominated by first derivatives, the preconditioning 
is carried out using central finite-differences on a refined mesh, as suggested in Section 5.2.2. 
The multi-dimensional preconditioner is inverted iteratively, using an ADI scheme related to 
the well-known Beam and Warming algorithm; fourth-order artificial viscosity is used in the 
preconditioner to remove odd-even decoupling of the central differences, and the ADI scheme 
uses pentadiagonal solvers to include the artificial viscosity terms. Early tests indicate that the 
preconditioned iteration has a spectral radius of about .2 on reasonable meshes, and that the 
overall time-stepping scheme is stable to Courant numbers of several hundred. 

8 Other Applications 

Spectral methods have been developed and applied to many other problems in viscous flow. 
This chapter will summarize some of the recent developments. 

CHQZ, Sec. 6.3 describe the use of spectral methods for solving the incompressible boundary- 
layer equations. This technique has been extended to the compressible boundary-layer equations, 
and has been used in the work of Macaraeg, Streett, and Hussaini (1988) in furnishing highly 
accurate self-similar mean flows for flat plate boundary layers and for free shear layers, both in 
supersonic flow. 

The linear stability equations of incompressible flow have been a very fruitful application for 
spectral methods (CHQZ, Sec. 6.4). Some interesting algorithmic issues are raised by the Orr- 
Sommerfeld equation, for it is a fourth-order equation and requires two boundary conditions at 
each endpoint. Bernardi and Maday (1988) have analyzed the alternative collocation schemes 
for fourth-order problems and have deduced which are the optimal formulations. They have 
studied domain decomposition versions as well. Phillips and Karageorghis (1989) present some 
results for a spectral domain decomposition method for a two-dimensional fourth-order problem. 
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A number of researchers have developed spectral codes for solving the compressible linear 
stability equations. Just within our group at NASA Langley, different spectral codes have been 
developed by Ng (1989), Macaraeg, Streett, and Hussaini (1988), Malik(1988), and Erlebacher, 
Herbert, and Hussaini (in preparation). Moreover, A. Thunn of the University of Stuttgart (1988, 
private communication) has developed a similar spectral code. These have been applied to jet, 
boundary-layer, free shear-layer, and plane Couette compressible stability problems. Macaraeg, 
Streett, and Hussaini (1988) utilized a staggered mesh in their algorithm. This avoids the need 
for a pressure boundary condition. Velocities and temperature were located on the Gauss- 
Lobatto points where the momentum and energy equations were collocated, and the pressure 
was located on the Gauss points where the continuity equation was enforced. A multi-domain 
version of the scheme, using the integral flux balance condition, was also implemented. It proved 
extremely useful in improving the efficiency of the spectral method, particularly in resolving the 
highly oscillatory modes that exist in supersonic flow. They presented comparisons with a 
second-order finite-difference scheme, Malik (1988) compared with a fourth-order method, and 
Thunn has made comparisons with up to eighth-order finite-difference methods. 

Spectral methods continue to be quite useful for investigations of secondary instability 
(CHQZ, Sec. 6.4). The reviews by Bayly, Orszag, and Herbert (1988) and Herbert (1988) 
provide a recent status report. Ng (1989) has used spectral methods to study the secondary 
instability of subsonic jet flows. 

Problems in reacting flow continue to challenge numerical methods. Guillard and Peyret 
(1986) describe a spectral method for a premixed flame problem in which a prescribed mapping 
is used to increase the resolution of the flame front. Bayliss, Gottlieb, Matkowsky, and Minkoff 
(1987) have pursued adaptive grid strategies for spectral methods. Augenbaum (1988) has 
studied adaptive solutions of the one-dimensional wave equation and the Burgers equation. 
Macaraeg and Streett (1988) have used a domain decomposition method to great advantage in 
studying non-equilibrium chemistry effects at hypersonic speeds. 

Several serious investigations have been made on the suitability of spectral methods for local 
memory parallel processors. These are, of course, more challenging for spectral methods than 
shared memory parallel processors due to the global nature of the approximation. Erlebacher, 
Bokhari, and Hussaini (1989) implemented on the FLEX/32 a three-dimensional explicit spec- 
tral code for simulating transition on a flat plate in supersonic flow. They achieved over 95% 
efficiency with as many as 16 processors. Fischer, Ho, Karniadakis, Ronquist, and Patera (1988) 
have run incompressible spectral-element simulations on a 16-processor Intel iPSC/2 hypercube 
at over 75% efficiency. Pelz (1989) has achieved over 80% efficiency for a 128 s simulation 
of isotropic turbulence using a pseudospectral method on a 1024-processor NCUBE. Finally, 
Tomboulian, Streett, and Macaraeg (1988) have implemented a three-dimensional spectral col- 
location algorithm for channel flow on a 16,384-processor Connection Machine 2, and achieved 
performance comparable to that obtained for the same algorithm on a single-processor of a 
Cray 2. Of all these applications only Pelz (1989) used FFT’s in his implementation. The 
others used matrix-multiplies to accomplish the Fourier/Chebyshev transforms and derivatives. 
The key to most of these implementations was the ability to perform fast matrix transposes 
and/or Fast Fourier Transforms in a manner that kept the calculation CPU-bound rather than 
communication-bound. 


50 



9 Trends 


A cursory analysis of the post-1985 citations in the bibliography indicates where current research 
on the development and analysis of spectral methods has been concentrated. Certainly the work 
on spectral domain decomposition methods stands out. Perhaps it has been a matter of necessity, 
but there is little question that the most innovative and demanding work on domain decom- 
position techniques in general is now being performed for spectral methods. The bibliography 
also suggests that there has been a resurgence of contributions to iterative methods, particularly 
spectral multigrid methods. Taken together these two trends portend an even broader range of 
application for spectral methods in the 1990’s. Moreover, recent achievements in implementing 
spectral methods on massively parallel computers indicate that spectral methods will keep pace 
with the on-going architectural revolution in computer hardware. 

We have made relatively few comments in these notes on the great strides that have been 
made during the 1980’s in the numerical analysis of spectral methods. The state-of-the-art has 
progressed from the ability at the start of this decade to analyze rigorously little more than linear 
model problems to the imminent capability of furnishing rigorous error estimates for spectral 
approximations to the full unsteady Navier-Stokes equations with general boundary conditions. 
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| Figure 1. Comparison of finite-difference (left) and Legendre spectral (right) differentiation, 

j The solid curves represent the exact solution, and the dashed curves are their numerical approx- 

imations. The solid lines are the exact tangents at x = 0, and the dashed lines the approximate 
tangents. The error in slope is noted, as is the number of intervals N. 
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Figure 2: Spectral solution of Taylor-Couette flow. The time history and frequency spectrum 
of axial velocity data is taken at 98% of the gap between the cylinders and on the mid-height 
plane. 
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Figure 3. Three sine waves which have the same k = -2 interpretation on an eight point grid. 
The nodal values are denoted by the filled circles. The actual sine waves are denoted by the 
solid curves. Both the k = 6 and the k = — 10 waves are misrepresented as a k = -2 wave 
(dashed curves) on the coarse grid. 
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Figure 4. Finite-difference (circles) and Fourier spectral (diamonds) approximations after one 
period to a periodic wave equations whose exact solution is represented by the curve. 
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Figure 5. Fourier collocation and second-order finite-difference solutions to the periodic Burgers 
equation problem at t = x/8. 



Figure 6. A non-overlapping domain decomposition for a grooved channel. 
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Figure 7. A overlapping domain decomposition for a grooved channel. 
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Figure 8. A one-dimensional domain decomposed into two subdomains. 
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Figure 11. Patched and overset grids for a two-dimensional hyperbolic problem. Heavy lines 
denote subdomain boundaries and light lines indicate the Chebyshev grids with the subdomains. 
(Courtesy of D. A. Kopriva.) 
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Figure 12. Convergence of spectral domain decomposition method for a hyperbolic system on 
the grid of Figure 11. (Courtesy of D. A. Kopriva.) 
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Figure 13. Eigenvalue spectrum of the Chebyshev first-derivative operator, preconditioned with 
central finite-differences on the spectral mesh, and on meshes refined by 20% and by 50%. 
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Figure 14. The multigrid V-cycle. The number of relaxations after restriction is denoted by Nd 
and the number of relaxations before prolongation is denoted by N u . 
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Figure 15. Convergence of spectral shock-fitting solution to the blunt-body problem. 
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Figure 16. Rotation (R), alternating (A), and Galerkin (G) results for turbulence intensity in 
homogeneous turbulence on 64 s , 96 s , and 128 s grids. 
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Figure 17. Alternative velocity and pressure nodes for flows with two non-periodic directions. 
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Figure 18. Contours of disturbance vorticity in the primary domain for the incompressible 
channel flow simulation. The solution in the buffer domain is not shown. 
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Figure 19. Surface disturbance vorticity on the lower channel wall for the incompressible channel 
flow simulation. The solid line is the simulation result and the dashed line is the linear theory 
prediction. 
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